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We discuss a phase space representation of quantum dynamics of systems with many de- 
grees of freedom. This representation is based on a perturbative expansion in quantum 
fluctuations around one of the classical limits. We explicitly analyze expansions around 
three such limits: (i) corpuscular or Newtonian limit in the coordinate-momentum repre- 
sentation, (ii) wave or Gross-Pitaevskii limit for interacting bosons in the coherent state 
representation, and (iii) Bloch limit for the spin systems. We discuss both the semiclassi- 
cal (truncated Wigner) approximation and further quantum corrections appearing in the 
form of either stochastic quantum jumps along the classical trajectories or the nonlinear 
response to such jumps. We also discuss how quantum jumps naturally emerge in the 
analysis of non-equal time correlation functions. This representation of quantum dynam- 
ics is closely related to the phase space methods based on the Wigner- Weyl quantization 
and to the Keldysh technique. We show how such concepts as the Wigner function, Weyl 
symbol, Moyal product, Bopp operators, and others automatically emerge from the Feyn- 
mann's path integral representation of the evolution in the Heisenberg representation. We 
illustrate the applicability of this expansion with various examples mostly in the context of 
cold atom systems including sine-Gordon model, one- and two-dimensional Bose Hubbard 
model, Dicke model and others. 
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I. INTRODUCTION 

Representing quantum mechanics using entirely the language of phase space variables, like 
in classical physics, attracted significant theoretical attention since its early days. It was clear 
from the beginning that such a description should be closely connected to statistical physics, 
where quantum fluctuations should be captured, at least partially, by statistical fluctuations. A 
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very 



important development in addressing this question came from the work o ' 



Moya l (jMoval 



,1949|) , who showed that the Liouville equation for the Wigner function (jWignerl . Il932l ) (Wigner 
transform of the density matrix) assumes the form identical to the Liouville equation for the density 
matr ix in the class i cal lim it if one generalizes the notion of Poisson brackets to Moyal brackets (see 



Ref. ( Hillery et al 



19841 ) for review) . This formalism closely relies on the correspondence b etween 



quan tum operators and classical functions of phase space variables first suggested by Weyl (jWeyll . 



19271 ): the expectation value of any operato r is equal to t 



operator weighted with the Wigner function (Hil lery et al. 



le av erage of the Weyl symbol of this 



19841 ). The Wigner function just plays 



the role of the probability distribution of phase space variables. Because the Wigner function is 
not positively defined, it is often referred to as the Wigner quasi-probability distribution. 

Phase s pace methods were a l so developed and found many ap plications in quantum optics (see 



e.g. Refs. (j Gardiner and Zollerl . 



2004; 



Walls and Milburn 



19941 )). There it is convenie nt to work 



1926) and then 



Glauber 



19631 ). 



with the quantum coherent states first introduced by Schrodinger (jSchrodinger . 
reintroduced into the quantum optics and termed as "coherent states" by Glauber 
Quantum classical correspondence is then found by associating creation and annihilation operators 
of photons with complex amplitudes, which also define the phase space variables. The Liouville 
equation for the coherent state Wigner transform of the density matrix takes a form very similar to 
the equation for the Wigner function in the coordinate momentum represent ation. More recently 



phase spac e 



Steel et al 



met hods in quantum optics were adopted to atomic systems (jSinatra et al. 



2002; 



1998) and a lready found numerous applications in the field of cold atoms (see e.g. 



Ref. (IBlakie et al. 



20081 ) for review). 



The Wigner function and the corresponding equations of motion give a natural framework for 
the representation of quantum dynamics using entirely the language of phase space variables. The 
time evolution in this language can be then represented through deterministic motion of phase 
space degrees of freedom accompanied by the quantum jumps. The main goal of this work is to 
discuss this representation in detail. Within this approach one naturally recovers the classical 
limit and the structure of quantum corrections. There are many classical limits we are familiar 
with from the elementary physics, each providing a natural choice of the phase space variables: i) 
corpuscular, particle, or Newtonian limit, where coordinates and momenta of particles constitute 
the phase space; ii) wave limit, where number and phase of the fields (or equivalently complex 
amplitudes corresponding to creation and annihilation operators) form the phase space. This limit 
naturally appears in optics, wave mechanics, or more recently in the physics of weakly interacting 
degenerate Bose gases; iii) mixed limit where some degrees of freedom are treated as particles and 
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some as waves. This limit naturally appears in gauge theories, most notably in electromagnetism; 
iv) classical spin limit, where spins are treated as rotators, and so on. The feature which unifies 
all classical descriptions is that the time evolution of the corresponding phase space variables is 
defined by the unique solution of deterministic equations of motion for given initial conditions. 
Very often even in classical systems the initial conditions are different from one experiment to 
another. Then one has to perform an additional averaging over the probability distribution of 
the initial conditions. Because the number of classical equations of motion scales linearly with the 
number of degrees of freedom, it is usually possible to simulate classical dynamics in relatively large 
systems. Contrary exact quantum simulations are usually limited to very small systems because 
of exponential dependence of the Hilbert space size on the number of degrees of freedom. 

In this work we will be interested in approaching non-equilibrium dynamics in quantum systems, 
where we start from a well defined initial state either pure or mixed, either equilibrium or not. Then 
we change the parameters of the Hamiltonian in a certain way, and follow the time evolution of the 
observables of interest describing the system. We assume that our system is closed (i.e. there is no 
external heat bath) and that there is a well defined Hamiltonian, possibly time dependent, which 
determines the time evolution. One of the very important questions, which needs to be addressed, 
is how the exponential complexity of quantum dynamics kicks in if we gradually increase h. More 
specifically one needs to address the following questions: i) what we should do with equations of 
motion, ii) what we should do with quantum observables, and iii) what we should do with initial 
conditions as we increase h starting from zero. We would like to note that the Planck's constant h 
plays the role of the saddle-point parameter such that in the limit S->0a unique classical theory 
is recovered. In practice, there can be other saddle point parameters which play the same role like 
the (inverse) mode occupation in bosonic systems, the (inverse) spin size in the spin systems and 
so on. So by h we effectively understand some parameter which governs the strength of quantum 
fluctuations. As we go along it will become clear what plays the role of this parameter near 
particular classical limits. 

These questions were partially addressed before. In particular, it was realized that the Wigner 



transform of the density matrix is especial 



when the quantum fluctuations are small ([Steel et al 



V suitable for ana l yzing dynamics near the c lassical limit 



1998; 



Walls and Milburn 



19941 ). It was also 



understood that in the leading order quantum fluctuations appear only through the (Wigner) 
distribution of the initial conditions while they do not affect the equations of motion themselves 1 . 



1 It should be noted that this statement is only true if the quantum Hamiltonian is written in the symmetric Weyl 
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The resulting semiclassical or truncated Wigner approximation (TWA) was 



various non-equilibrium problems with bosonic cold atom systems (see Ref . (jBlakie et all . l2008l ) for 



successfully 



applie d to 



review). There were also developed schemes based on other representations of the density matrix 



in th e coherent state basis, in particular, based on the positive P-representatiq n (jDrummond et al. 



2007 



Gardiner and Zollerl . 



2004 



Steel et al 



1998 



Walls and Milburn 



19941 ). Somewhat related 



method allowed one to map proble m of interacting lattice boso ns into sto chastic evolut ion using 
the basis of Fock (number) states (jCarusotto and Castinl . 120031 ) . In Ref. (jDeuarl . |2009| ) a mixed 
representation of the density matrix, which smoothly interpolates between positive P- and Wigner 
representations, was suggested. Such schemes, however, have their o wn issues and do not generally 



allow for systematic expansion in quantum fluctuations. In Refs. (jPolkovnikov 



2003a 



bj) it was 



shown that TWA naturally appears in the coherent state path integral representation in bosonic 



systems as the semiclassical approximation incorporating quantum flu ctuation in t 



r e leadin g order 



to the classical (Gross-Pitaevskii) dynamics. Furthermore in Ref. (jPolkovnikov 



2003al ) it was 



shown that the consequent quantum corrections to the expectation value of a particular observable 
appear in the form of the nonlinear response of this observable to the infinitesimal changes in 
the classical fields across the trajectories. Recently TWA was also extended to o ther situations 



going beyond traditional domain of interacting spinless Bose systems. E.g. in Ref. (jAltland et al 



20091 ) the TWA was used to analyze 



In Ref. (jBerg et al. 



functions. In Ref. ( 



the Dic ke model. In Ref . (ISau et al 

Ml 



dynam ics of a particular spin-boson system described by 



2009) TWA was first applied to the spinor condensates. 



2009) for the first time TWA was extended to tackle multi-time correlation 



Mathey and Polkovnikovl . 



20091 ) TWA was applied to analyze non-equilibrium 



vortex formation due to quantum fluctuations in a two-dimensional quantum rotor model. 

The same semiclassical approximation (though not termed as TWA) was developed in the 
context of the quantum dynamics in t he co o rdinate-mom entum representation using the language 



of the Moyal brackets (jHillery et al 



1984 



Zurek 



2003). In the leading order in h, where the 



Moyal brackets reduce to the Poisson brackets, one finds an analogue of TWA, where quantum 
fluctuations enter only through the initial conditions while the equations of motion are affected 
only in higher orders in H. This approach was exte nsively used f or studying quantum chaos and 



decoherence in single-particle systems (see e.g. Refs. (|Habib et al 



Zurek 



1998; 



Karkuszewski et al 



2002; 



2003)), particularly for a single particle system coupled to a heat bath. It was argued that 



form. Otherwise there can be non-vanishing corrections to the equations of motion, which remain, however, fully 
deterministic. 



6 



the decoherence in chaotic systems was responsible for extending the validity of the semiclassical 
approach by eliminating effects of interference between single particle trajectories. In the single- 
particle context TWA was also applied to analyze collisions of wave packets with one-dimensional 



barriers ( Lee and Scully 



19831 ). 



In driven many-particle systems the most com monly used approach to the non-equilibrium phe- 



nomena i s the Keldysh diagrammatic technique (jKamenevl . 



Keldvsh 



2002 



Kamenev and Levchenko 



2009; 



19651 ). This technique is often used to derive the quantum kinetic equation and ob 



tain single-particle distribution functions in non-equilibrium systems (jKamenev and Levchenko! . 



200J). 



Usually this technique is applied for finding steady sta t es in driven systems (see E.g. 



Refs. (Dalla Torre et al. 



, 2009; 



Mitra et al. 



2006; 



Takei and Kim 



20081 )). b ut it can be also used 



for a n alyzing dynamics in situations, where initial conditions are important (jAltland and Gurarid . 



2008 



Altland et all |2009j). While formally the Keldysh technique and phase space methods are 
equivalent, there is an important difference between them 2 . The elementary objects in the phase 
space methods are the Wigner function and phase space variables. In the Keldysh technique the 
elementary objects are the single particle's Green's functions. Closely related to the Keldysh di- 
agrammatic technique there are also non-equilibriu m approaches based on the functional integral 



representation of evolution (|I.Plimak et al. 



2001a 



Plimak et al. 



accompanied w i th various approximat i ons like effective actions 



Dias and Pratal . 



2007; 



Gasenzer 



2009; 



Rev et al 



2001|) , which can be further 



Braunss 



2010 



Dias et al. 



2006; 



, 120051 ) . These methods are now actively devel- 



oped. 

In this paper we will discuss how quantum dynamics can be formulated entirely in the phase 
space using the language of deterministic (classical) trajectories and stochastic quantum jumps. 
We will review both the earlier results and present the new ones. A very important goal of this 
work is to show that the formalism is basically identical in various phase space representations 
(coordinate- momentum, coherent state, and spin coherent state). While some details between 
these representations are different the structures of the description are very similar suggesting 
generality of this approach. We will review how all these results can be derived using the Feynman 



path integral formulation o f the time evolu t ion. T 



in the Keldysh technique ([Kamenevl . 



2002, 



2005 



lis derivation resembles the one routinely used 



Kamenev and Levchenko 



2009J) with the only 



difference that instead of working in the Schrodinger representation, where the density matrix 



2 There is some ambiguity in defining what the Keldysh techniques is. Sometimes this techniques is understood in 
the broad er sense of the functio nal representation of the time evolution on the Schwinger-Keldysh contour. See 
e.g. Ref. jl.Plimak et al.ll2001al ') 
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evolves in time, we will work in the Heisenberg representation, where the operators describing 
observables change in time. In the coherent state representation the methods d escribed in this 



work are also related to the functional integral approach used by Plimak et al. (ll.Plimak et al 



200 lal ). There the density matrix only determines the initial state. The main difference of that 
work with the approach discussed here is that we will not attempt to find the evolution of the 
whole density matrix, which is an exponentially complex object in many-particle systems and 
contains tremendous amount of unmeasurable information. Instead we will be interested in finding 
expectation values of particular observables, like the energy or its fluctuations, various equal or 
nonequal time correlation functions, etc. We will show that such concepts like the Wigner function, 
Weyl ordering, TWA, and quantum corrections naturally appear from the Feynmann's path integral 
without need to make any assumptions. We will show that quantum corrections to TWA can be 
written either in the form of the nonlinear response to infinitesimal quantum jumps or in the 
form of stochastic quantum jumps with non-positive quasi-probability weight. In both cases each 
jump carries a factor of the effective Planck's constant squared. The representation of the quantum 
operators through the phase space va riables also naturally emerges fro m the path integral formalism 



in the form first suggested by Bopp (jBopp 
dependent operators: 



1961 



Hill erv et al 



x(t) -> x(t) + 



ih d 
~2 dp(t) ' 



or 



x(t) -)• x(t) 



ih d 



pit) -> p(i) 



p(t) -> pit) + 



1984) and generalized here to time 
ih d 



2 dx(t) 



ih d 



(1) 



(2) 



2 dp(ty w ^ w 2 dx(ty 
where the partial derivatives are understood as the response to the infinitesimal quantum jumps 
occurring at time t and the choice of the sign is dictated by casuality of the evolution (see the 
next section for details). This Bopp representation will play the key role in the whole formalism. 
Throughout this paper we will reserve "hat" -notations to denote quantum operators and "non- 
hat" -notations for the functions defined in the phase space . There is a representation similar to 
Eqs. (H]) and ([2]) for the creation and annihilation operators (see Eqs. (|72p and (|73p ). Without 
further details it is clear that the representation ([I]) of the quantum coordinate and momentum 
operators above immediately reproduces the classical limit at H — )■ 0. It is interesting to note 
that the requirements of the casuality of the representation of the quantum dynamics through the 
phase space trajectories dictates that only certain multi-time correlation functions are allowed (see 
Sec. E]). These tu rn out to be precisely the same correlation function s which appear in the theory 



of measurements (Clerk et al 



2008; 



Nazarov and Kindermann , 



2003) 
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Let us note that in three (and higher) dimension al continuous systems with short range in- 



terac tions TWA suffers from ultra-violet divergencies (IBlakie et al 



2008 



Deuar and Drummond 



20071 ). This indicates that the expansion of dynamics around the classical limit can be ill defined 
and likely requires some short-distance renormalization of the bare parameters of the Hamiltonian. 
In the context of cold atoms certain ways to deal w ith these divergenci es by effectively truncating 



20081 ). However, comprehensive 



the Hilbert space were suggested in the literature (jBlakie et al. 
understanding of these issues is far from being complete and is beyond the scope of this work. 
It is not even understood whether these ultra-violet divergencies are intrinsic property of three- 
dimensional systems or they are simply related to the singular na ture of short-range poten tials 



in three dimensions, which is well known from equilibrium physics ([Fetter and Waleckal . Il97ll ). A 
simple way to avoid these divergencies, as we illustrate in Sec. WN\ is to work with lattice systems. 

The paper is organized as follows. First in Sec. [U we give a brief introduction into the phase 
space representation of the density matrix and operators. Starting from the coordinate-momentum 
representation we describe such concepts as the Wigner function, Weyl symbol, Moyal product and 
Moyal brackets, and Bopp operators. We then discuss similar concepts in the coherent state picture. 
In particular, we introduce coherent state analogues of the Poisson and Moyal brackets as well as 
of the Moyal product. In the end of this section we discuss close analogy between the coordinate- 
momentum and coherent state representations. Then in Sec. IHII we give a detailed overview of the 
formalism discussed in this paper. We tried to highlight all the features of the formalism important 
for practical implementation of TWA and quantum corrections to particular problems skipping the 
details of the derivations of these results. We split this section into three parts discussing first 
dynamics in the coordinate- momentum representation, then in the coherent state representation, 
and finally quantum dynamics of spin systems. Then in Sec. IIVI we discuss the implementation of 
the formalism to particular problems. We start from the simplest possible example of a harmonic 
oscillator suddenly driven from equilibrium gradually increasing the complexity of the systems and 
describing dynamics in the sine-Gordon model, one and two dimensional Bose-Hubbard models, 
Dicke model and others. For small systems we give comparison of the results obtained within this 
approach with the exact results and discuss in which situations the quantum expansion works well 
and in which it fails. Our goal was not to describe the corresponding physics of these models in 
detail, but rather to focus on the applicability of the method. Then in Sec. |V]we guide the reader 
through the derivation of the results using the path integral approach. This section is rather 
technical and can be skipped unless the reader is interested in understanding these details. One 
of the main purposes of this section is to show that the Weyl symbol, Wigner function, quantum 



9 



corrections, casual representation of the dynamics, Bopp operators, etc. naturally appear from 
the path integral approach without need to make any assumptions or to justify this choice a- 
posteriori. Finally in Sec. I VII we discuss connections of the path integral method used here and the 
other approaches based on the von Neumann's equation for the density matrix. In Appendix [Al 
we illustrate the detailed implementation of the first quantum correction to the TWA for the 
sine-Gordon and Bose-Hubbard models. In Appendix [B] we discuss a somewhat subtle issue of 
emergence of the Weyl ordering of the Hamiltonian in classical equations of motion within the 
path integral representation of the evo l ution (which usually requires the normal ordering of the 



Hamiltonian (IKamenev and Levchenkd . 



2009|)). 



II. BRIEF OVERVIEW OF THE PHASE SPACE REPRESENTATION OF QUANTUM SYSTEMS. 



In this section we will briefly introduce the basic tools, which will be essential in the next sections 
when we will talk about dynamics. In particular, here we will introduce such concepts as the Wigner 
function, Weyl symbol, Moyal bracket, and Bopp operators. The Moyal brackets will also allow us 
to write down equations of motion for the density matrix in the phase space representation and 
discuss their close connection to the classical Liouville's equations (see Sec. IVip . We will discuss 
and contrast the two basic phase space representations: coordinate-momentum and the bosonic 
coherent state. The first representation is naturally connected to the corpuscular classical limit 
while the second representation corresponds to the wave classical limit. In this work we will not 
talk about the Grassmann's number representation of the fermionic phase space. 



There are several reviews availab 



quantum systems. In particular Ref. (Hi llery et al. 



e in literature desc ribing the phase space representation of 



19841 ) gives an excellent overview of such rep- 



resentation in the coordinate-momentum phase s pace. Detai l s of the coherent state represen- 



tation of bosonic system s can be found in Refs. (jBlakie et al 



Walls and Milburn 



2008; 



Gardiner and Zoller 



2004; 



19941 ). We will thus not attempt to give a comprehensive review of these rep- 
resentations and skip the proofs, which can be found in the above mentioned texts. We will only 
introduce the tools necessary for understanding the consequent sections describing dynamics. We 
will also contrast these two phase space representations. By introducing coherent state Poisson and 
Moyal brackets we will show that these two descriptions of quantum systems identically map to 
each other under the change of the phase space variables. The reader familiar with these concepts 
can skip this section and directly proceed to the next one. 



Let us note that the phase space representation of quantum systems is not 



10 



unique (j Gardiner and Zollerl . 12004k iHillerv et all 1 19841 ; IWalls and Milburnl . 1 19941 ). In this work 



we will focus only on the Wigner-Weyl representation because, as we will show later, it automati- 
cally emerges from the path integral description of quantum dynamics and gives the most natural 
connection between classical and quantum dynamics. This representation also treats conjugate 
phase space operators symmetrically. 

A. Coordinate-momentum representation 



In the context of ordinary canonical variables like coor dinates and momenta the phase space 



approach to quantum mechanics was pioneered by Moyal (jMoval 



1949 ). It i s based on the Wey l 



ordering of operators an d the Wigner distribution function (see e.g. Refs. (Hi llerv et al. 



Wang and Connell . 



1984: 



198£j)). 

Weyl symbol and Wigner function. The key ingredient in the phase space description of quantum 
systems is the so called Weyl symbol, which gives a one to one map between quantum operators 
and the ordinary functions defined in the phase space. For Hermitian operators this map is real. 
Thus for an arbitrary operator Q(x,p) the Weyl symbol f2jy(x, p) is formally defined as 



x + 2 / ex P 



(3) 



We use the vector notations to highlight that we are dealing with general d-dimensional multi- 
particle phase space of the dimension 2D, where the factor of two reflects that for each degree of 
freedom we are dealing with pairs of conjugate variables. In particular, for N-particle system in 
d-dimensions we have D = Nd. If the operator (l(x, p) is written in the symmetrized form then it 
is straightforward to see that the Weyl symbol £lw is obtained by simple substitution x — > x and 
p — > p. In particular, this is true for all operators of the form 0(x, p) = -A(x) + B(p). 

If the ordering in f2(x, p) is such that the coordinate operators appear on the left of momentum 
operators then the Weyl symbol ([3j) can be also written as 



(4) 



Here f2(x, p) is the function obtained from the operator (1 by direct substitution x — > x and 
p — > p. The equivalence of Eqs. (|3|) and (jl]) can be established by inserting the identity resolution 
i = 1/2/ dr} |p + f//2)(p + 77/2| into Eq. 

The second key ingredient of phase space methods is the Wigner function which is defined as 
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the Weyl symbol of the density matrix p: 



W(x,p) 



P 



x+|^e ip ^ ri = J d£p(x-£/2,x + £/2)e ip ^ r \ 



For any proper density matrix the Wigner function is normalized such that 

dxdp 



r W(x,p) 



1. 



(5) 



(6) 



{27Th) D 

The Wigner function together with Weyl symbols of various operators gives complete description 
of a given system. In particular, the expectation value of any operator Q(x, p) is found by averaging 
the Weyl symbol of this operator over the phase space with the Wigner function playing the role 
of the probability distribution: 

dxdp 



<n(x,p)> 



;W(X,P)JV(X,P)- 



(7) 



{2irh) D 

Because the Wigner function is not positively defined it is often referred to as the Wigner quasi 
probability distribution. 

M oyal product . W eyl symbols of operators satisfy important Moyal product rela- 



tion ( Hillerv et al. 



19841 ) which defines the Weyl symbol of the product of two operators 



through the Weyl symbols of the individual operators: 
(fiin 2 )w(x,p) = fii,w-(x,p)exp 

where 



ih . 
—A 
2 



A 



d d do 
dp cbc cbc dp 



fi2,w(x,p) 



dpj dxj dxj dpj 



(8) 



(9) 



is the so called simplectic operator. It is straightforward to verify the validity of Eq. (jSJ) using 
the definition of the Weyl symbol ([4"]). For any two functions A(x,p) and B(x, p) the product 
— AAB = {A, B} defines the Poisson bracket between these functions: 



{A B} - V 9A 9B - 9A 8B - {A B} 
^-r 1 dxj dpj dpj dxj 



(10) 



Note there are different sign conve ntions in d e fining Poisson brackets; e.g. compare 



refs. (Landau and Lifshitz 



of Rcf. (Hill erv et al 



1982 ) and (jfflllerv et al 



1984 ). In Eq. (JTUJ) we adopted the choice 



19841 ) . So the Moyal product of the two operators (J8|) can be also expressed 



through the exponent of the Poisson bracket: 

(^i0 2 )tv(x,p) = Oi )W (x,p)exp 



ih 



{...} 



^2,w(x,p). 



(11) 
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From Eq. ([8]) it follows that the Weyl symbol of the commutator of the two operators ti, 
can be expressed as 



n 



A 



where {. . . } MB defines the Moyal bracket: 

{• • • }mb 



^2 w = 2i£l\ w sin 



{•••} 



2W 



= ih{ni,Q 2 } M B > 



f^X) ^2 



(12) 



sm 



{...} 



(13) 



In the limit ft. — > the Weyl symbol of the commutator of any two operators reduces to the Poisson 



bracket between the corresponding classical function s multiplied by ih 



known manifestation of the correspondence principle (jLandau and Lifshitz 



ris is of course a well 



1981 



)• 



Bopp operators. Wigner-Weyl quantization is intrinsically co nnected with the symmet ric rep 



resentation of the phase space operators first suggested by Bopp (|Boppl . 



1961 



Kubol . 



196J): 



x 



ih d 



P = P 



ih d 
Tdx"' 



(14) 



It is clear that this representation satisfies the canonical commutation relations. There is an 
equivalent Bopp representation based on the left derivatives: 



x = x 



ih *d 
2~dp' 



P = P + 



ih *d 
2~dx" 



(15) 



where the left arrow implies that the derivative acts on the operator on the left. Loosely speaking 



the representation (|15|) is obtained from Eq. (|14p by integrating by parts. The forma! 



that the Bopp operators re produce the Weyl symbol can be found in Refs. (jHillerv et al. 



proof 



Osborn and Molzahn 



1984; 



2002 ) . We will also prove these statement in Sec. [V] when we will analyze 



non-equal time correlation functions. 

Let us illustrate how the Bopp operators can be used in practice. Suppose we are interested 
in a Weyl symbol of the operator xp. For simplicity we choose a single particle in one dimension. 
Then using the Bopp operators we can immediately construct the Weyl symbol of this product: 



xp = xp 1 



x + 



ih d 



P 



ih d 
~2dx~ 



1 = xp + 



ih 



(16) 



2 dp J V" 2 dx J r ' 2 
Thus we conclude that {xp)w = xp + ih/2. The same result can be obtained using the right 



derivatives: 



xp 




ih 

xp + —. 



(17) 
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We inserted unity into Eqs. (|16p and (|17p to illustrate that the are no other operators on the right 
of p and left of x. The same result can be obtained directly from Eq. as well as from the 
symmetrization 

_ xp + px 1 ih 

xp= h -[x,p\ -> xp+ — (18) 



and from the Moyal product expansion: 

ih 



(xp)w = xexp 



A 
2 



p = xp — —xAp = xp + — . (19) 



which is exactly equal to the result obtained from Eq. (j4|). As we will see below in the Heisenberg 
picture the derivative terms in Eqs. (|14p and (|15j) get an additional interpretation as a response to 
infinitesimal quantum jumps at a particular moment of time. Then the choice of the left or right 
derivative is dictated by the casuality of the evolution. 

B. Coherent state representation. 

All the results above immediately generalize to the coherent state representation of the phase 
space. Bosonic coherent states \ip) c are defined as eigenstates of bosonic annihilation operators t/j 
such that t/)\ip} c = ip\i>)c- Explicitly the coherent states read: 

|V) c = exp[^V t ]|0) = V^=|n), (20) 

where |0) is the vacuum state and \n) = {^) n /^/n\ |0) is the n-particle Fock state. It is straight- 
forward to check that 

^)c = ^)c (21) 

Obviously there is a direct analogy between the action of the coordinate x and momentum p 
operators on the coordinate basis (x\x) = x\x), and p\x) = —ihd x \x)) and the action of the 
annihilation ip and creation ^ operators on the coherent state basis. As we will see this analogy 
will persist throughout the paper. The coherent states (120p are not orthogonal and not normalized. 
In particular 

C (V#') C = exp^V]- (22) 

Because of non-orthogonality, the coherent state basis is over-complete. Using these states one can 
write the resolution of the identity 

1= f dil>dil>*e-M*\il>) c M, (23) 
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where the integration measure is defined as dipdtp* = ^HipdQip/ir. 

Weyl symbol and Wigner function. Using coherent states one can define the Weyl symbol of 
and arbitrary operator Q^,^') as 



tt w {ip,%p+ 



2 M 



drfdi) ( tp 



tp + 1 > (■ 



(24) 



Here the vector ip = {ipj} consists of complex amplitudes corresponding to different single-particle 
eigenstates, M is the Hilbert space size. So j can denote coordinates, momenta, different internal 
spin states, etc. As in the coordinate-momentum representation the Weyl symbol of f2 is the partial 
Fourier transform of the matrix elements of this operator between different coherent states. The 
additional gaussian factors come from the normalization of the coherent states. Again as in the 
coordinate momentum representation the Weyl symbol of a symmetrically ordered operator can be 
obtained by simple substitution ■?/>—)• tp and — > ip*. For normally ordered operators Eq. ([24 
gives 



(25) 



Similarly to the coordinate-momentum representation the Wigner function is defined as the 
Weyl symbol of the density matrix: 



W{ip,xp* 



1 

w 



drfdr] ( ip 



P 



?7 

1p + - )r 



(26) 



The expectation value of any operator, by analogy with Eq. ([7|), is given by the average of the 
Weyl symbol with the weight given by the Wigner function: 



dipdip* w(ip, *p*)n w (ip, v*)- 



(27) 



So the Wigner function again plays the role of the quasi-probability distribution of the complex 
amplitudes. 

Moyal product and Bopp operators. The Weyl symbol of the product of two operators can be 
written by analogy to Eq. ([8]): 

~A r: 



(flii} 2 )w = ^i,w exp 



2W, 



(28) 



where we introduced the symplectic coherent state operator (cf. Eq. ([9])): 

d o d a 



a c = E 



dipj dip* dip* dipj 



(29) 
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Similarly to the coordinate-momentum representation it is convenient to introduce the coherent 

state Poisson bracket between two arbitrary functions A and B as 

9 A dB dA dA 
{*,B} c = AA c B = J:——-——. (30) 

j 3 3 J 

The coherent state Poisson bracket is a classical analogue of the commutator. In particular, 
{ipi,ipj} c = is a classical analogue of = 8ij. In general it is easy to see that in the 

classical limit of large occupation numbers [A, B] — > {A, B} c . Note that the coherent state Poisson 
bracket defined above (Eq. (|30p ) directly corresponds to the commutator without need to multiply 
by ih. This choice is dictated by the convenience so that the classical Hamiltonian equations of 
motion in the wave limit read ihdtipj = {ipj,H} c . It is easy to check that for bosonic systems 
with two-body interactions this equation becomes identical to the Gross-Pitaevskii equation (see 
Sec. [VI] for details). 

Using the coherent state Moyal product rule ([8]) and the symplectic structure of the operator A c 
we can find the expression for the Weyl symbol of the commutator of the two operators Q = [Oi, SI2] 



&2,W = ^2,w} MBC i (31) 



= 2^1^ sinh 

where we introduced the coherent state Moyal bracket by analogy with standard Moyal bracket (|12p : 

1 



{■■■)mbc = 2sinh 



„<■■■>. 



(32) 



Wigner Weyl quantization is intrinsically associated with the coherent state Bopp operators 
which can be introduced by analogy with Eqs. (Q2 



(33) 

* i ^* i+ *wr* i ~*w (34) 

The choice of the representation with the conventional (right) derivatives and the one with left 
derivatives is arbitrary. As we will see below for time dependent problems it is dictated by casuality. 
This representation of creation and annihilation operators is clearly symmetric and preserves the 
correct commutation relations. It also automatically reproduces the Weyl symbol of any operator 
(the proof of this statement will be given in Sec. IV.Bj) . To illustrate this let us show a couple of 
examples of computing Weyl symbols of anc ! ^ ^ ipip first using the Moyal product rule and 
then using the Bopp operators. Thus using Eq. (|28|) we find 




(^) w = ^(i_±_^_^_ )^ = |V> 2 |- ~. (35) 
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The same expression can be obtained using the Bopp operators 

W = «i - (** - \£) (* + \^) 1 = W - 5" <36> 

Note that the Weyl symbol for the number operator can be obtained also by symmetrizing it first 
ip^ip = l/2|$T^j + ^>t] — 1/2 and then substituting t/> — > if) and ^ — > tp* in the symmetrized form. 
This Weyl symbol can also be obtained by direct application of Eq. (|25p to the operator Cl = "ip'y). 
Let us also briefly consider the other example first using the Moyal product: 

mm)w _ {r? {, _ I^g + J**) ^ _ |, T _ + I (37) 
and then using Bopp operators 

*H^)V^) 2l ^ 2|2 - 2| * |2 4 ( 38 > 

Again both methods give the same result which can be also verified using the direct integration in 
Eq. ([25]). 



^ ^ 'pip 



C. Coordinate-momentum versus coherent state representations. 

To summarize the discussion above we will contrast the two phase space pictures in Table HI 
The main purpose of doing this is to emphasize the similarity between the coordinate-momentum 
and coherent state representations. Normally the latter is introduced in the context of the sec- 
ond quantization. This table shows that there is no secondary aspect to the quantization in the 
coherent state representation. These two representations correspond to two dual descriptions of 
the phase space: corpuscular and wave. Note that this dual description is only available for 
bosonic s ystems. For fermion i c coh erent state phase space description requires use of Grassmann 



variables (jNegele and Qrland 



1988), which do not have a natural classical interpretation. While 
one can formally define the Wigner-Weyl quantization for the fermionic systems, its extension to 
dynamics is not understood by now and will be a subject of the future research. 



III. QUANTUM DYNAMICS IN PHASE SPACE. 

In this section we will discuss the representation of quantum dynamics in the phase space pro- 
viding sufficient details for further exploration of this formalism and for its practical applications. 
We will skip the details of the derivation of these results, which will be discussed later in Sec. El 
First in Sec. IIII.Al we will overview the situation in the coordinate-momentum representation, where 
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TABLE I Coherent state versus coordinate momentum phase space 



Representation 


coordinate-momentum 


coherent 


Phase space variables 


x,p 




Quantum operators 


x,p 




btandard representation 


x —> x, p > —ihd x 
(coordinate basis) 


(coherent state basis) 


Canonical 
commutation relations 


(a, t3 refer to different particles) 


[V'i,^]] = <*y 
(i, j refer to single-particle states) 


Quantum-classical 
correspondence 


x ->• x, p ->• p, [A, B] ->• B} 


/ 4 Rl — V 9A dB dA dB 
X^i^ic j dipj dip* dip* dip j 


tit „ .. r i • 

vVigner tunction 




p 


x+ &\ e ip£/H 


WW, *!>*) = | //d»i*dt| (V-!|/5|V + !) 


Weyl symbol 


fiw(x,p) = Jd£ (x- f 


n 


x+ 0e* pe / fi 


✓ /, i ■ \ 1 f* f* i j, j / j Til If . TJ \ 

-0 ) = UI dr i*dv (V> - §1 n |V> + f) 

x e -M 2 -~lrM 2 e^'?*^-^*) 


Moyal product 


(r2ifi 2 )vK = ^i,iyexp [-^A] fi 2 ,iy, 
A — V ^ ^ 

^-^a dp a dx a dx a dp a 


(fiifi 2 )w = Oi^exp [4f] n 2) w, 

V 1? V_ if 

c — <9v>* &i>i 9ipj 


Moyal bracket 


{ni,n 2 } MB = -|Oisin[|A] n 2 


{Q 1 ,Cl 2 } MBC = 2fii sinh [±A C ] ft 2 


Bopp operators 


» _ j_ a a _ ih *d 
P = P-#J; = P + #i 


^ = </>*" 5^ =</>*+||;, 



the phase space is represented by the coordinates and momenta of particles. In the corresponding 
classical limit the particles move according to deterministic trajectories in this phase space de- 
fined by the microscopic Newton's (or more generally Hamiltonian's) equations of motion. Next in 
Sec. IIII.Bl we will describe dynamics of interacting bosons in the coherent state phase space. There 
the classical (Gross-Pitaevskii) limit corresponds to the (matter) waves. And finally in Sec. IIII.CI 
we will discuss spin systems. Using the Schwinger mapping of spins to bosons and then using 
the results of Sec. IIII.BI we will avoid the need of working with rather complicated spin-coherent 
states. As we already mentioned in this paper we will not talk about the wave limit for fermions. 
In many situations dynamics of fermions is described by bosonic collective excitations and thus 
can be analyzed within the methods described here. We will mention an example of such mapping 
of fermions to bosons and the application of the phase space formalism to the dynamics when we 
discuss the Dicke model in Sec. IIV.GI 

Before proceeding let us note that there is a quite common misconception that the Gross- 
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Pitaevskii equations are only applicable to the condensates and describe only the mean-field dy- 
namics, i.e. that they apply only to situations when the condensate is described by a single wave. 
This statement is equivalent to saying that the Maxwell's equations in electromagnetism are only 
applicable to lasers (coherent states of photons) or that the Newton's equations are only applicable 
to rigid bodies when the motion of all particles can be described by a single center of mass degree 
of freedom. In reality the Gross-Pitaevskii equations are simply classical wave equations of motion 
similar to the Newton's equations, Maxwell's equations, or Bloch equations. This point perhaps be- 
comes clearer as we go along and see t hat the Gross-Pitaevskii dy namics can be expressed through 



the standard Hamiltonian formalism (ILandau and Lifshita . 



1982l ) using the language of coherent 



state Poisson brackets introduced earlier in Sec. HH 
A. Coordinate-momentum representation 

1. Semiclassical dynamics of a single particle in an external potential. Truncated Wigner Approximation. 

Let us start the discussion from considering a single particle in one dimension moving in an 
external, generally time dependent, potential V(x,t). In quantum mechanics the time evolution of 
a wave function (density matrix) is described by the Schrodinger (von Neumann) equation: 

ih?^^=H(t)*(x,t), (39) 



or 

_dp(x, x' , t) 



ih- 



H(t),p(x,x',t) , (40) 



dt 

where ^(x,t) (p(x, x',t)) stand for the wave function (density matrix) of the particle; square 
brackets denote the commutator and 

H(t) = ^ + V(x,t) (41) 

is the Hamiltonian. We use the capital \& for the wave function to avoid confusion with coherent 
state complex amplitudes. The corresponding classical Hamiltonian equations of motion are given 
by the standard formulae: 

^ = {*,«}, % = {PW, (42) 

where 



{A, B} = d x Ad p B - d p Ad x B 
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is the Poisson bracket. The classical equations of motion should be supplemented by the initial 
conditions. Together they determine a unique classical trajectory. However, it is important to em- 
phasize that in principle there can be some statistical uncertainty in knowing the initial conditions. 
This is especially true in many-particle systems. In this case one needs to average the observable 
of interest over these initial conditions. For example, the average of some classical observable 
Q c \(x,p, t) at some moment of time t is given by 

(Q c l(x,p,t)) = J dx dp W c i(x ,p )n c i(x(t),p(t),t), (43) 

where W c \(xo,po) is the probability distribution of the initial coordinate and momentum. In the 
fully deterministic case xq = x and po = p we have W c \(xo, po) = 5(xo — x)S(po — p) and we are 
back to the unique value of O with no averaging needed. 

In quantum mechanics the coordinate and the momentum of a particle can not be simultaneously 
defined due to the uncertainty principle. Therefore it is not possible to fully localize the particle in 
the phase space and thus the distribution function should always have nonzero width. A natural 
object which substitutes the classical probability distribution £l c \(x,p) in the quantum case is the 
Wigner function or equivalently the Wigner transform of the initial density matrix introduced in 
Sec. [TI] (see Eq. ([5])). In the purely classical limit h — > the Wigner function reduces to the classical 
probability distribution. In particular, for the equilibrium density matrix at finite temperature the 
classical limit of the Wigner function is the Boltzmann's distribution. Similarly the object which 
replaces the classical observable is the Weyl symbol of the quantum operator (see Eq. Q). 

Interestingly in the leading order in quantum fluctuations Eq. ()43|) . which expresses the expec- 
tation value of any operator as an average of the Weyl symbol weighted with the Wigner function, 
immediately generalizes to the time dependent problems. Namely, up to the second order in h 
quantum fluctuations do not affect the classical equations of motion, which now play the role of 
the characteristics along which the Wigner function is conserved (see also Sec. IVIj) . This is equiva- 
lent to the classical picture, where the probability distribution W c \(x(t),p(t),t) is conserved along 
the classical trajectories. Thus in the leading order in quantum fluctuations Eq. (|43p still remains 
valid if we replace W c i(xo>Po) ~~ ► W(xo,po) and Q c i(x(t),p(t),t) — > Qw(x(t),p(t),t): 

(Cl(x,p,t)) « J J dx dpoWo(xo,po)Q w (xci(t),p c i(t),t), (44) 

Comparison of Eqs. (|43|) and (|44p highlights a close connection between the quantum and the clas- 
sical statistical averaging. Since, as we noted earlier, the Wigner function is generally non-positive 
its direct interpretation as a probability distribution becomes somewhat problematic. For equilib- 
rium initial ensembles at high temperatures the Wigner function coincides with the Boltzmann's 
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function and thus becomes positive. However, as the temperature is lowered the difference between 
classical and quantum distributions becomes stronger: the classical distribution gets narrower both 
in terms of coordinates and momenta while the quantum Wigner distribution always satisfies the 
minimum uncertainty principle. The semiclassical approximation (|44p to the dynamics is very 
important because it gives the first step in going from classical to quantum description of the 
dynamics. A similar approximation in th e coherent state phase spa ce (see below) was termed as 



truncated Wigner approximation (TWA) (jWalls and Milburnl . 1 1994! ) for the reasons which will be 



clear in Sec. IVII We will stick to this terminology also in the coordinate-momentum representation, 
thought it is not commonly used there, because this is essentially the same approximation. 

2. Nonequal time correlation functions. 

The semiclassical (TWA) approximation can be extended to finding non-equal time correlation 
functions like 

(n 1 (x,p,ti)n 2 (x,p,t 2 )). (45) 

In the classical limit such correlations do not depend on the ordering of the operators. However, 
quan tum mechanically the ordering is important. Within the phase space approach (see also 



Ref. (jBerg et al. 



20091 )) this non-commutativity can be elegantly absorbed into the language of 
quantum jumps, which naturally emerge extending the notion of Bopp operators (j!4j) and f)15|) to 
the Heisenberg representation: 

iW ^ w + l4 = lW "i4 (46) 

« f ^<<>-2^=^> + 2&Ki)' (47) 
Now the derivatives above are understood as quantum jumps, i.e. if we are interested in measuring 
e.g. (x(ti)x(t 2 )) with t\ < t 2 then at the moment t\ the conjugate momentum undergoes an in- 
finitesimal jump: p{t\) p(ti) + 5p{t\) . After that the system continues to evolve deterministically 
(within TWA) and at the moment t 2 in addition to the expected term x(ti)x(t 2 ) there is an extra 
contribution equal to the response of x(t 2 ) to this jump 5p{t\) multiplied by ih/2. I.e. 

{x(h)x(t 2 )) « J J dxodp Wo(xo,po)x(t 1 )x(t 2 ) + j J J dx dp W (x ,p ) > ( 48 ) 

where in the last term the limit 5p{t\) — > is implied. The expression above is approximate only 
to the extent of TWA being approximate. If t\ > t 2 then one can still use the same procedure 
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introducing the jump at t\. But then in order to evaluate the response of xit^j it becomes necessary 
to propagate the equations of motion backwards in time. This is inconvenient and unphysical, since 
it violates the casuality of the description of the dynamics. Instead one can use the left derivative 
rule, introduce the jump in pfa) and evaluate the response of x(t\) with respect to this jump 
multiplied by —ih/2: 



(£(ii)x(t 2 )) 



dx dp Wo(x ,p )x(ti)x(t2) 



ih 



dxQdp Wo(xo,po 



5x(tx) 

] s P (t 2 ) 



(49) 



In this way we clearly preserve the casuality of the evolution. It is interesting to note the higher 
order correlation functions including more than two different times can be comp uted in a casual 



way o nly for special orderings of the operators (see also the discussion in Ref. (ll.Plimak et al 



2001a 



)). We will return to this point in Sec. (jV]) when we discuss the non-equal time correlation 
functions in detail. 

From Eqs. (|46|) and (|47p it is obvious that if we are interested in finding the expectation value 
of the symmetric combination of say x(ti) and xfo) or x{t\) and p{t2) then the quantum jumps 
cancel each other and we can use classical substitutions x — > x and p — > p to find such averages. 



For example, 



x(ti)x(t 2 ) + x(t2)x(ti) 



X(t!)x(t2). 



(50) 



Thi s is consistent wi th general expectations that symmetric correl ation functions are classi 



cal (IClerk et al 



2008 ) which are measured by ideal classical detectors (jNazarov and Kindermann 



2003). In general for multi-time averages purely classical substituti ons occur if we ar e interested in 



time-symmetric ordered correlation functions (see Sec. |V] and Ref. (jBerg et aZ.1 . 120091 )). Conversely 
the antisymmetric correlation function has only the contribution coming from quantum jumps: 





x{tl)x{t2) — x(t2)x(t\) — > ih; 



x(t 2 ). 



(51) 



dp{h) 

Clearly in this case extending Bopp operators to Heisenberg representation according to Eqs. 
and (|47|) also gives a natural extension of the commutation relations to non-equal time operators 
through the response to quantum jumps. In fact this statement is exact and is not limited to the 
TWA. 



3. Beyond the truncated Wigner approximation: nonlinear response and stochastic quantum jumps. 



The truncated Wigner approximation is formally exact if we are dealing with harmonic systems. 
For a single particle this means that the potential V{x) should be quadratic. If this is not the case 
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then TWA is only approximate. There are two equivalent ways of representing quantum corrections 
to TWA both using the language of quantum jumps. In the first way the corrections are represented 
as the nonlinear response (of the observable of interest) to these jumps and in the second way these 
jumps are stoch astically distributed. The s econd representatio n is so mewhat reminiscent to the 



Langevin noise ([Gardiner and Zollerl . 



2004; 



Walls and Milburn 



1994 ). however, with important 



differences. The representation of quantum dynamics through the noise also appears in quantum 



optics if one deals with two-body potentials and writes the Liouvi 



matrix in the P- or posit ive P-representation (jDrummond et al 



Walls and Milburn 



2007; 



l e equation for the density 



Gardiner and Zoller 



2004; 



19941 ). However, the jumps in the P-representation do not give a systematic 



expansion in quantum fluctuations. In the formalism discussed here based on the Wigner-Weyl 
quantization each jump carries an extra factor of ft? and thus the expansion in the number of jumps 
is equivalent to the expansion in the powers of the Planck's constant. As we show next the jumps 
have several unusual properties and only superficially resemble standard noise. 

For simplicity we will consider again a single-time expectation value of some observable Cl. 
This simplification is not really important. In a more general case one can combine Eqs. (146 p and 
(|47p with the results presented below to compute quantum corrections to multi-time correlation 
functions. Rather than giving a complete general formula we will show the first few terms in the 
quantum expansion: 



(n{x,p,t)) 
i - 



dx Q dpoWo(po,xo) 



* 1 ft 2 d 3 V{x) d 3 

o T 3!22 1? dx(r) 3 9p(r) 3 
t rt 



* , 1 h*d 5 V(x) d 5 
dr- 



o 5!2 4 i A dx(r) 5 dp{rf 
ft 4 d 3 V{x) d 3 d 3 V(x) d 3 



+ L ^ L ^ (3! 2 2 ) 2 * 4 &c(ri) 3 dp{ Tl ) 3 dx(r 2 ) 3 dp(r 2 



+ ... \n w (x(t),p(t),t). (52) 



The expansion above is clearly in powers of ft 2 . In the leading (zeroth) order one simply recovers 
TWA. In the next order in ft? one needs to introduce a quantum jump in the momentum: p(r) — > 
p(r) + 8p(r), calculate the third order nonlinear response of the observable to this jump, multiply 
this response by the third derivative of the potential with respect to x(r) as well as by other factors 
and integrate over r. In the next order of ft 4 one can either have two jumps with the third order 
response or a single jump of the fifth order response. In general each jump of 2n + 1-th order 
carries a factor of ft 2n . All the factors in Eq. ([52]) are written in the way easily generalizable to 
the higher order terms. In Sec. IIVI we will show some examples illustrating how finding these 
quantum corrections can be implemented in practice. From Eq. (I52p it is evident that TWA is 
asymptotically exact at short times. The smaller ft the longer the regime of validity of TWA. 
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At long times quantum corrections may or may not become important depending on the details 
of the problem and the observable f2. Moreover the moment where TWA breaks down can have 
very different sensitivity to the actual value of h. We will see this considering various examples in 
Sec. IIV1 Thus it is very hard to make a general statement of the time of validity of TWA based 
on Eq. (I52p without knowing details about the system. 

There is another formal representation of Eq. (1521) through the stochastic quantum jumps: 



(n(x,jp,i))« / J dx Q dpoWo(x ,Po) 



+ 



(2 2 ) 2 



... \n w (x(t), P (t),t). 



(53) 



YlEhihlil_ Similarly to the ordinary Langevin 



To simplify notations we introduced V2 n +i,i — — Q x ( T .ji n +i 
diffusion we discretize the time here in steps of the size At: Ti = iAt. The notation Spi = ^vAr 
above implies that at the moment Tj the momentum undergoes a quantum jump of the magnitude 
£i \/ At where the £j 6 (—00,00). The weight function Fs(£i) can be interpreted as a (quasi)- 
probability distribution of this jump. With this interpretation the jumps become stochastic. We see 
that in the leading order in quantum fluctuations at some moment in time during the evolution the 
momentum p(r^) undergoes a single stochastic jump. The distribution function of the magnitude 
of this jump i*3(£) must satisfy the requirements that all its moments are finite and the moments 
up to the second vanish: 

/oo roo roo /*oo 

F 3 (£K = 0, / £F 3 (Od£ = 0, / fF 3 (OdH = 0, / £ 3 F 3 (£R = 1. (54) 
-00 J— oo J —00 J -co 

Similarly in the order ffi we encounter either two third order jumps or a single fifth order jump. In 
the latter case £ is distributed according to the function i^(£) where first four moments are equal 
to zero. In general i 7 2n+i(£) must satisfy 

/OO /'OO 
r\F 2 n+i(£R = 0, m<2n, / f n+l F 2n+1 (^ 
-00 J — OO 

The functions i^n+iXO are not unique. One possible choice is 



1. 



(55) 



F, 



2n+l 



(0 



1 1 - r JL] e -e/2 



2«+i v ^F(2n + l)! if2 " +1 y^2 



where H2 n +i is the Hermite polynomial. In particular 

1 /£ 3 A 1 



-e/2 



(56) 



(57) 
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The representation (|53p bears some similarities with the diffusion but there are important 
differences too. Thus in ordinary diffusion jumps are always proportional to V At, while in Eq. (|53|) 
the power of the time interval At is smaller. Likewise in usual diffusion only the first moment 
of the distribution function should vanish, which allows one to use the Gaussian function for the 
-F(£). In Eq. (|53H at least the first two moments of F(£) vanish so this function can not be positive 
definite and can not be interpreted as a classical probability of the jump. As in the nonlinear 
response representation of quantum corrections, each jump explicitly carries a factor of ft? or 
higher (depending on the order of jump) allowing for a consistent expansion of time evolution in 
quantum fluctuations. We note that while the representation (|53p is formally equivalent to the 
one given by Eq. (|52p in the limit At — > 0, practical convergence of these two expansions can be 
quite different. In particular, as we will show using several examples, the step At can be chosen 
to be much bigger than the one determined from the requirement that the two trajectories with 
and without jump remain close to each other throughout the full evolution. This is also true 
about usual diffusion where the distribution function is much more robust to the choice of the 
time step than individual trajectories. Let us mention that for the coherent state pha se space the 



repre s entation of the ev o 



1998 



Gevorkvan et al 



integration ( 



ution through the cu b ic nois e was introduced ear 



1997; 



I.Plimak et al 



Gardiner and Zollerl . 



200J). 



2001al : 



Plimak et al. 



ier in Refs. ( Gevorkvan 



200ll ) using Ito's stochastic 



4. Generalization to many particles and higher dimensions. 

The formalism above can be straightforwardly generalized to interacting many-particle systems 
in arbitrary dimensions. Coordinates and momenta become vectors and acquire an additional 
particle number index. For example Eq. (|52p generalizes to 



(Cl(x,p,t)) ai J f dx dp W (po,x ) 



ft i fi? d 3 V(x) d 3 \ 

Here indices a, f3, 7 run over different components of the D-dimensional phase space, e.g. for a 
system of iV particles in three dimensions a, /?, 7 = x±, y±, Zi,... xn,Vn, zn- As usually we imply 
summation over repeated indices. This expression can be again rewritten through the quantum 
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diffusion. Thus instead of Eq. (|53|) we get 



(n(x,p,t)) ~ J J dx Q dp W (x ,p ) 

1 + f^/ /n^ £ af we |« d ( a! (t),p(t),*),(59) 

L i 3 3 m <r(a,P, 7 ) a P 7 ^(r^^^/A^J 

where a"(a, /3, 7) stands for all non-equivalent permutations of indices a, f3, 7, e.g. a < (3 < 7. The 
functions F a ^ n (£) should satisfy the conditions that all its moments smaller than two vanish and 

f fl[dUUp^Fa,pM = l' (60) 
3 3 m 

A possible choice for the functions F a p~ is 



Fa,a,a ) 



2^2^ V 3 



4tt 



(61) 
(62) 
(63) 



In Eqs. (|61|) - (|63|) a, /3, and 7 stand for non-equal indices and no summation over identical indices 
is implied. 



B. Coherent state representation. 

As we pointed in the previous section for bosonic systems, or more generally the systems which 
support collective bosonic excitations, it can be more convenient to expand around the wave classi- 
cal limit and to use the coherent state description of the phase space. Examples of classical waves 
include phonons, photons, Bogoliubov excitations in superfluids, matter waves etc. In this section 
we are going to address the issue of representing quantum dynamics in this phase space. 

Rather than keeping a completely general discussion we will concentrate on bosonic systems 
with two-body interactions. In particular, such systems appear in the context of cold atoms. 
The analysis itself, however, does not rely on this assumption and can be applied to arbitrary 
n-body interactions. As we will see below the structure of the quantum dynamics in the coherent 
state phase space is identical to that in the coordinate momentum phase space. Let us point 
that the emerging sem iclassical approximation or TWA was first adopted to interacting bosons in 



Ref. (Steel et al 



19981 ). Since then there were many works applying the TWA to p articular cold 



atom experiments. For the review of some of these applications we refer to Ref. (jBlakie et all 
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2008). The corrections to TW A in the form of nonlinear response were first obtained by the author 



in Ref. (jPolkovnikovl . 



2003a|). 



Let us consider the following Hamiltonian: 

U = J2 h i$Hi + \ E » ( 64 ) 

j ijkl 

where iM and ifj are the bosonic creation and annihilation operators. The indices i,j,k,l can 
refer to any single-particle basis, e.g. they can describe coordinate or momentum states or the 
single-particle states in an external potential. The Hamiltonian is split into the noninteracting 
part with the matrix elements h®j and the interacting part characterized by the matrix elements 
of the interaction potential Uijki- Both h° and u are allowed to depend on time. The operators ^ 
and if) satisfy the canonical commutation relations: 

$ i J]]=8 ij . (65) 

The Hamiltonian (|64j) is quite general. For example, if indices i,j denote spatial positions on the 
lattice by the appropriate choice of the matrix elements of h° and u it reduces to the Hubbard 
model: 

W-hub = - J ^(Mi>j + h.c.) + E v jnj + 2 E U v : ™ ifl i : ' ( 66 ) 

3 V 

where n, = ipjipi is the density operator, Vj is the external potential, J is the tunneling matrix 
elements between nearest neighbors, and Uij is the matrix element of the density-density inter- 
action. We use capital letters to distinguish parameters of the Hubbard model from those in the 
general Hamiltonian. The semicolons in the Hubbard Hamiltonian imply the normal-ordered form 
of the interaction, i.e. the form where the creation operators appear on the left of the annihilation 
operators. We note that in the continuum limit, when the lattice spacing goes to zero, the Hub- 
bard model reduces to that of interacting bosons without the lattice. The indices ... in the 
Hamiltonian (j64j) can also refer to momentum states. Then in the absence of external potential 
q = €k5 pq , where is the energy of a single p article with momentum k and Vk+q :P -q tPj k are the 



Fourier components of the interaction potential (jPethick and Smith! 



2001 



As we discussed in Sec. |H] the classical wave limit in the coherent state picture corresponds to 
large occupation numbers of different modes. In this case we can treat the operators ip^ and ijj 
as complex c-numbers. In the classical limit one substitutes the commutation relations (I65p with 
the coherent state Poisson brackets (i30j) : {ipi,ip*} c = 5ij. Under this substitution the Heisenberg 
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equations of motion for the fields tpj, ijA become the Gross-Pitaevskii equations: 



ih 



dt 



&Hw 

dip* 



^2 + ^2 Uijkq^j^klpq, 

3 j,k,q 



(67) 



where 



^ij — h%j a ^ ] ^ikkj V'ikjk ^kijk ^kikj ■ 



(68) 



The difference between hP and h comes from the Weyl ordering of the Hamiltonian. 

For example for the Hubbard model (see Eg. (1661)) the equation above reduces to the fa miliar 



discrete Gross-Pitaevskii equation (jPolkovnikov et al. 



2002; 



Trombettoni and Smerzi 



2001 



ih 



dt 



(69) 



where 



Vi 



(70) 



Here Oj denotes the nearest neighbors of the site j. We see that for the Hubbard model with 
translationally invariant interactions (Uij = Ui+k,j+k for arbitrary k) the Weyl ordering only gives 
a term proportional to the total number of particles, which can be absorbed into the overall phase 
or the chemical potential. However, as we will illustrate in Appendix [Bl the extra terms can be 
very important in other situations with e.g. spatially dependent interactions. Let us note that the 
Planck's constant in Eq. (I69j) can be safely set to unity by either redefining energy or time units. 
This is in fact custo mary in cold atom systems where relevant energy scales are often measured 



in Hz (see e.g. Ref. (jBloch et all [2008)). The real parameter, which determines the degree of 
quantum fluctuations is the occupancy of relevant bosonic modes: classical waves correspond to 
high occupation numbers. 



1. Truncated Wigner Approximation 

Similarly to the coordinate-momentum picture in the leading semiclassical order (or TWA) 
quantum corrections appear only through the initial conditions and the form of the observable 

M,ft,t))= f J d^d^w (^Mn w (^(t),^*(t),t), (7i) 

The integral above is taken over different realizations of all components of initial values of if) and 
if}* so that dtpodtp* stands for Ylj dipojdip*j. The fields tp(t) and tp*(t) are related to tpo and 
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i/)q via the solution of Eqs. (|69|) . The Wq{x]jq, i/)q) is the Wigner transform of the initial density 
matrix ([2BJ) : and Ow(i/>(f), i/)*(t), t) is the Weyl symbol of the operator 0, §M$ . 

For computing non-equal time correlation functions we can use the recipe similar to the one 
discussed above in Sec. IIII.AT21 In particular, we can extend the coherent state Bopp operators (j33|) 
and (I34p . to the Heisenberg representation: 

%®^®-hm-w t)+ *<>m- (72) 

* w ^* w+ i^ = * w -bjfw (73) 

interpreting the partial derivatives as response to infinitesimal quantum jumps. Namely, the deriva- 
tive with respect to tpj(t) in the equation for V>j(i) implies that at moment t the field ipj(t) undergoes 
an infinitesimal jump ipj(t) —> ipj(t) + 5ipj(t). Then one needs to evaluate the derivative of all op- 
erators appearing on the right of ^\(t) with respect to this jump. The derivatives above can also 
rewritten in terms of real and imaginary parts of ifi according to the standard rules: 

d Id id d Id id 

~d~4> ~ 2~imp ~ 2~dlhp : d^* ~ 2£m^ + 2 5^0' ^ ' 

2. Quantum corrections. 

As in the case of the coordinate-momentum representation, quantum corrections to Eq. (|7ip 
can be written either in the form of nonlinear response or stochastic quantum jumps. In the first 
way (cf. Eq. ([58]) ) we get 

n w (^(t),i>*(t),t). (75) 




ijkl 



d 3 

^ {T) d^(T)d^ k (T)d^(T) ~ C ' C - 



The interpretation of this expression is also similar to what we encountered earlier. In the leading 
order in quantum fluctuations beyond TWA one has a single quantum jump where the classical field 
i^j (r) undergoes the infinitesimal transformation ipj (r) — > ipj (r) + 8ipj (r) . Then this field continues 
to evolve according to the classical equations of motion. In the end one calculates the nonlinear 
response of the desired observable to this jump. Further corrections appear as multiple quantum 
jumps. Note that for the two-body interactions there are no higher order quantum jumps which 
would correspond to fifth and higher order derivatives with respect to tp and ip* (cf. Eq. (j58[) ). As 
we discussed earlier, near the classical limit fields ip have large occupation so the derivatives with 
respect to tfj and ip* give factors of the inverse square root of the occupation number. 
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One can represent these corrections also using stochastic quantum jumps: 



i-^EE^ / dT %i 



n ijkl 



n w (ij(t),^(t),t) 



(76) 



where a = j,k,l and dT^ kl is the phase space volume element with jkl denoting all nonequivalent 
permutations of indices i,j,k. For j ^ k ^ I this phase space volume element is given by dT^ kl = 
d£jd£* d^kd^dfydQ ; forj^k = l it is dT^ kk = d£jd£jd£kd£k and so on. The interpretation of the 
jump 5ip a in Eq. (|76p is the same as in the coordinate-momentum representation: at the moment 
T n the complex field ip a (T n ) jumps to ip a ( T n) + 6ifi a (T n ). The function F can be interpreted as the 
quasi-probability distribution of the stochastic quantum jump. It should satisfy the requirement 
that all its moments up to the second vanish and the third moments are equal to: 



k-i ) 



(77) 



As in the coordinate-momentum case these requirements do not define F uniquely. A possible 
choice for this function is 



^(0,6,6) = #&6e-tel a -K*l a -tol 3 , j^k^l, 
^',£;,&) = My 2 -l) e" 
F(^,C k ,Ck) = C-f k e-^ 2 -^\ 3 + k 
n^Mi) = £(l£/-2) e 



For the Bose-Hubbard model with the Hamiltonian (1661) and local interactions U, 



(78) 
(79) 
(80) 
(81) 



U5ij the 



general formula (|76j) reduces to 



iU 



/ j d^l^ir^F^j) - c.c.)n d (m,^(t),t) 



(82) 



where as in the Eq. (|76|) at the moment r n the classical field ipj undergoes a quantum jump 
ipj -)• 4>j + £j VAt and F(£,-) is given by Eq. (|8T|) . 



C. Spin systems. 



In this section we will discuss expansion of dynamics for interacting spin systems or more 
generally dynamics of systems with angular momenta. We will work in units where h = 1. The 
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classical limit corresponds to S, size of the spin, being large: S 3> 1 so the parameter of the 
expansion around the classical limit should be 1/5. The spin operators satisfy the canonical 
commutation relations: 



(83) 



where e a b c is the fully antisymmetric tensor. Let us note that fo rmally spin dyn amics can be 



mapped to the dynamics of bosons using Schwinger representation (jAuerbachl . 



1998 ): 



a^a 



crp, s = fi'a. 



(84) 



This representation allows us to apply results from the previo us section directly to the spin s ystems 
without need to introduce spin-coherent states (see e.g. Ref. (jTakahashi and Shibatal . Il975l )). The 
bosonic fields a and (3 in Eq. f|84j) should satisfy an additional constraint n = cVa + j3* j3 = 2S. 
In the imaginary time formalism of the equilibrium statistical mechanics this constraint poses an 
important problem and should be taken care of introducin g an auxiliary constraint field at each 



moment of time (see e.g. Ref. (jArovas and Auerbachl . 12009 )). However, in real-time dynamics the 



situation is much less complicated because this constraint is preserved in time. Indeed h com- 
mutes with all spin operators and thus with any Hamiltonian with arbitrary spin-spin interactions. 
Therefore if the constraint is satisfied at the initial time then the dynamics of spins is equivalent 
to the unconstrained dynamics of bosons under the mapping (|84p . Thus semiclassical approxima- 
tion (TWA), quantum-classical correspondence, and quantum corrections for spins can be directly 
deduced from the results of the previous section. 

Quantum classical correspondence Using Eqs. (J72J) and (j73|) we can find an analogue of the Bopp 
operators for the spin systems: 

a* a - p*/3 If d 2 d 2 



•5 + 



8 \ da* da 
d_ 
da 




13- 



<9/3*<9/3 
1 d 2 



a 



d 
da* 



a 



d_ 

da 



4 dad (3* ' 
1 d 2 



4 da*d/3 ' 

These equations can be also written using compact notations: 



% 

s 

2 





1 


s x V 


~ 8 







V + (s • V)V - -sV 2 



where V = d/ds. For example, for s z this correspondence gives 



d_ 

ds„ 



d 

' ds x 



id s z f d 2 _ dP_ 

8ds~ z ~ 16 yds 2 ~ ds 2 



ds 2 



d 2 



d 2 



8 ds T d8, 



8 ds y ds 2 



(85) 
(86) 
(87) 

(88) 
(89) 
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These formulae can be used in constructing Weyl symbols for various spin operators. Let us 
give a few specific examples: 



1 ' 



(90) 



Note that the presence of i in the latter expression is not surprising since the product s z s x is not a 
Hermitian operator. Like before with coherent states the spin Bopp operators can be also used to 
find non-equal time correlation functions if we interpret derivatives as a response to infinitesimal 
jumps. 

Wigner transform. In principle, the mapping (|84p is sufficient to express the Wigner function 
of any initial state in terms of the bosonic fields a and (3. General expressions can be quite 
cumbersome, however, one can use a simple trick to find a Wigner transform of any pure single 
spin state and then generalize it to any given density matrix. Assume that a spin is pointing along 
the z-axis. This can always be achieved by a proper choice of a coordinate system. Then in terms 
of bosons a and (3 the initial state is just |25 + 1, 0), where S is the size of the spin. In other words 
the wave function is a product of two Fock st ates one having 2S + 1 par t icles and one part icles. 



The corresponding Wigner transform is then ([Gardiner and Zoller 

W(a, a*, (3,(3* 



2004 



Polkovnikovl . 



2003aT ): 



2e- 2 H 2 -2|/3| 2 L25+l(4 |a| 2 ) 



(91) 



where Ln(x) is the Laguerre's polynomial of order N. At large S the Laguerre polynomial is a 
rapidly os cillating function an d the Wigner transform is strongly localized near \a\ 2 = 25 + 1/2 
(see Ref. (jPolkovnikovl . l2003al )). So in this case to a very good accuracy (up to 1/S 2 ) one can use 



W(a,a*,/3,/3*) « v^e-^l^d 



a 



25-1/2). 



(92) 



At the same time the fluctuations in (3 are usually much more important since they represent 
fundamental uncertainty between S z and the transverse components of the spin. Reexpressing the 
Wigner transform f)91 1) in terms of spin components then we get 

2 



W(s z , s ± ) = —e-^L 2S+ Ms z + s]), 

7TS 



(93) 



where s± = \/s 2 + s 2 and s = \ s 2 + s 2 , represent transverse component and the total magnitude 



J x ' "y 

of the classical spin respectively. The Wigner function is normalized in such a way that 

f2w 



ds z \ sj_ds ± / d4>W(s x ,s±) = l, 

-oo J J 



(94) 
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where <f) is the polar angle of the spin in the xy plane. When S is large we have s z 
peaked around S and then Eq. (|93j) reduces to 

1 



W(s z ,s ± ) « — e 



- /S 5(s z -S). 



s strongly 



(95) 



This Wigner function has a transparent interpretation. If the quantum spin points along the z 
direction, because of the uncertainty principle, the transverse spin components still fluctuate due 
to zero point motion so that 



(4) = (si) = I 



(96) 



This is indeed the correct quantum- mechanical result. Clearly from Eq. (|93p one can derive the 
Wigner function for a spin with an arbitrary orientation by the appropriate rotation of the coor- 
dinate axes. 

For more complicated correlated spin systems the trick with aligning all spins along a particular 
axis does not work. In this case one can always rely on the Schwinger representation to get the 
Wigner function and then reexpress this function through the spin components . Alternatively one 



can directly work with the spin coherent states (jTakahashi and Shibata 



19751 ) and compute the 



Wigner transform without referring to the bosonic representatio n of spins. An example of such 



computation for two spin one half particles can be found in Ref. ([Cohen and Scully 



1983). 



Quantum corrections. Let us consider only the situation with two spin interactions of the form: 

where a,b = x, y, z. The corrections to the TWA can be again characterized by quantum jumps. 
These jumps can be recast in the form of nonlinear response so that within the first order in 1/5 
we have the following expression for the expectation value of an arbitrary observable: 

r * 



(Q(t)) = JdsoW(so) l+l JdTj^Jijhi^x V4(2(s i (r).V i )V 6J - S6J (r)V|) 



n w (s(t),t), 

(98) 



where Vj = d/dsj(r) and V&J = d/dsbjir). This expression should be understood in the same 
way as e.g. Eq. (|75p for the coherent state evolution: at the moment r the spins undergo in- 
finitesimal jumps sj(t) — > s, (t) + 5s j and at the end one evaluates the nonlinear response of the 
observable to these jumps. From the equation above it is clear that 1/S serves as the expansion 
parameter determining the strength of quantum fluctuations. Actually the small parameter is 1/S 2 
since the classical limit is achieved when 5 — > oo, J — > such that JS =const. 
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IV. APPLICATIONS 



Let us now discuss applications of the formalism described in the previous section to a number of 
sample problems starting from a single-particle moving in a harmonic potential and then gradually 
increasing complexity towards interacting many-particle systems. For implementing TWA and 
quantum corrections one essentially needs two ingredients: calculating the initial Wigner function 
and solving corresponding classical equations of motion. For complicated initial states finding 
the Wigner function can be a very challenging task by itself. Since in this paper we are mostly 
concerned with understanding various aspects of the dynamics we will put the discussion of finding 
initial Wigner function aside by starting from simple initial states, which correspond to equilibrium 
in systems with a quadratic Hamiltonian. In many situations complications associated with finding 
the initial Wigner function can be avoide d by either doing y arious approximations to the initial 



2008 ) or by starting from a simple initial 



state like the Bogoliubov's approximation (jBlakie et all 
state and then adiabatically evolving t he Ha miltonian in time so that the couplings reach the 
desired values ( Polkovnikov and Wang . 12004 ). This method is closely related to experimental 



procedures used in cold atoms, where one usually cools the system first and then adiabatically 



drive s the system towards the desired regime e.g. by turning on optica" 



2005 



2008) . This method is a l so closely related to quantum annealing algorith ms 



Farhi et al 



2001 



Kadowaki and Nishimori . 



1998; 



Santoro et al 



lattice 



Bloch et al 



Das and Chakrabartil . 



20021 ). If the evolution is 



sufficiently slow then one moves along the constant entropy curve. As long as one remains in the 



regime of sm all quantum fluctuations this method of adiabatical 



very reliable (jBistritzer and Altmanl . 
to note that in some low dimensional 



2007; 



Polkovnikov and Wand . 



y evo 



ving the Hamiltonian is 



20041 ). Though we would like 



slow rates (jPolkovnikov and Gritsevl . 



syst ems the adiabatic limit can be achieved at anomalously 



20081 ) . After the desired initial parameters of the Hamiltonian 



are reached via this procedure one can start real dynamical simulations. 



A. Single particle in a harmonic potential. 

As the first example let us consider a single particle moving in a harmonic potential. The 
advantage of starting from this simple situation is that all calculations can be done explicitly 
analytically. The Hamiltonian of a single harmonic oscillator is 

Ho = ^ + ^fx 2 = hw^H + 1/2), (99) 
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where the coordinate and momentum operators x and p are related to creation and annihilation 
operators ip and in a standard way: 



(100) 



mw y \ zn \ muj 

Now suppose that the particle is prepared in the ground state and we are suddenly applying a 
linear potential V(x) = —Ax. So that the Hamiltonian becomes 



% =%q — \x 



(101) 



Next we compute various observables as a function of time. 

Coordinate-momentum representation. First we will solve this problem using the coordinate- 
momentum representation. The ground state wave function for a harmonic oscillator is: 



ip (x) 



1 



(27r)V4^ 



-z 2 /4a2 



(102) 



where ao = y / h/2muj is the oscillator length ([Landau and Lifshita . Il98ll ) . According to Eq. ([5]) 
the Wigner transform corresponding to this wave function is 



W(x Q ,po) = J dS4>*(xo + Zm(xo - H/2)e ip °t /h = 2exp 



2al 



(103) 



where go = H/2oq = W mcoh/2. In this case W(xq,po) is positive definite and can be straightfor- 
wardly interpreted as the probability distribution of the initial coordinate and momentum. Next 
we need to solve the classical equations of motion: 



dp 
~dt 



-moj 2 x + A, 



dx p 
dt m 



satisfying the initial conditions x(0) = xq, p(0) = pq. Clearly the solution is 

x(t) = x c \(t) + xq cosiut) + —5- sin(wt), 

muj 



(104) 



(105) 



where x c \(t) = A/mw 2 (l — cos(w(t))) is a classical result describing the motion of the particle which 
is initially set to rest. Then we need to substitute this solution to the observable corresponding to 
the quantum operator of interest and find the average over the initial conditions. 

For the expectation value of the positi on we trivially find {x(t)) = x c i(t), which is just a 



particular case of the Ehrenfest's principle (jShankar 



19941 ). Similarly we find 



(x 2 ) = x 2 (t) = x 2 cl (t) + a 2 . 



(106) 



This is of course also the correct result, which can be easily obtained from the solution of the 
Schrodinger equation. 
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Next let us show how to compute a non-equal time correlation function. In particular, (x(t)x(t')) 
with t < t'. In this case according to Eq. (I46p we need to substitute x(t) — > x(t) + ih/2d/dpt and 
proceed with the classical calculation. This substitution is equivalent to an infinitesimal quantum 
jump in the momentum: p(t) — > p(t) + 5p. At the later moment t' we need to evaluate the response 
of the coordinate to this jump. Then we find 



(x(t)x(t')) = ( x c i(t) + xocos(ujt) + mxiut) + 
\ moj 



2 dSp 



(p dp 
Xd(t) + cos(cjt') H sin(wi') H sin(uj(t' — t 
moo muj 

= x c i(t)x c \(t') + Oq cos(uj(t — t')) + iog sm(u)(t' — t)). (107) 

Note that this correlation function is complex because it does not correspond to the expectation 
value of a Hermitian operator. Similarly for the correlator with the opposite ordering of t and tf 
we find 

{x(t')x(t)) = x cl (t)x cl (t') + a 2 cos(w(t - t')) - ia 2 sin(w(t' - £)) (108) 
There fore the symmetric part of the correlator corresponding to the result of a classical measure- 



ment (IClerk et al. 



20081 ) is given by 

/ X(t>)x(t) + X(t)x(t>) \ = + a 2 CQs{u;{t _ 0) m 

and the quantum antisymmetric part is 

(x(t)x(t') - x(t')x(t)) = 2ial sm{oj{t' - t)). (110) 

The quantum part vanishes in the limit H — > as it should. 

Coherent state representation. Now we will illustrate how the same results can be obtained in 
the coherent state picture. In the second quantized form the Hamiltonian of the system reads 

% = tiw($xf) + 1/2) - Aa o (-0 + (Ill) 

The classical equation of motion for the complex field ijj can be obtained from Eqs. (|30p and fjl 1 1 [> : 

dip 

ih-— = huip — Aao, (112) 
ot 



which has the solution 



l p(t) = ^ (1 - e~ iut ) + Voe"^. (113) 
nui 
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Using Eq. (|26p one can show that the Wigner transform of the initial vacuum state is 




(114) 



Now we can find the desired expectation values 



(x(t)) =aofoKt)+V*(*)) 



2a 2 \ 

TiIjJ 



(1 — cos(cjt)) = X c \(t). 



(115) 



Similarly 



(x 2 (t)) = og(v 2 (t) + mt)) 2 + www) = + 4 



(116) 



We obviously got the same answers as before. Similarly one can verify the result for the non-equal 
time correlation function. Of course it is not surprising that both methods give identical exact 
results for harmonic systems. However, it is important to realize that once we deal with more 
complicated interacting models the correct choice of the phase space can significantly simplify the 
problem. Moreover the expansions around the two possible classical limits are very different. Thus 
for a system of noninteracting particles moving in some external potential TWA in the coordinate- 
momentum representation is only approximate unless the potential is harmonic. At the same time 
TWA in the coherent state representation is exact. 

B. Single particle in a Mexican-hat potential. 

Let us now move to a somewhat more complicated problem, where we will still deal with a 
single-particle physics but in a quartic potential so that the Hamiltonian is: 



We assume that initially k(Q) = and the particle is in the ground state. Then n(t) changes in 
time according to 



At the moment when n(t) becomes negative classical equilibrium becomes unstable and the particle 
should move to one of the new minima (see the inset in Fig. [1]). The parameter 5 in Eq. (|118p 
controls the rate of the change of the Hamiltonian. Note that classically the particle with the initial 
conditions x = 0, p = will always remain at the same position. However, zero point fluctuations 
allow the particle to fall down to one of the minima and form the "Schrodinger cat" state. In 




(117) 



n(t) = < 



I -St 0<t<2/5 
-1 t>2/5 



(118) 
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FIG. 1 Time dependence of the variance of the displacement of a particle initially prepared in the ground 
state of the Hamiltonian (|117|) with k = 1 and A = 0. Then the potential evolves according to Eq. (|118|) with 
5 = 0.1. The other parameters are A = 0.1, m = 10, ujq = 1. We work with the units where ft = 1. The solid 
black curve represents the exact solution of the Schrodinger equation and the dashed red line represents the 
result of TWA. The inset shows the evolution of the shape of the potential with time. 

Fig. [T] we plot the variance of the coordinate as a function of time for the parameters 5 = 0.1, 
A = 0.1, m = 10, and u)q = 1. It is easier to work with the units, where H = 1. Then it is the 
inverse mass which plays the role of the quantum parameter: instead of h = 1 and m = 10 we 
could choose h = 0.1 and m = 1 and get identical results up to rescaling of units. As it is clear 
from the figure, the motion near x = indeed becomes unstable after n{t) becomes negative and 
(x 2 (t)} rapidly increases. Then at a certain point t = 2/5 the potential becomes stationary and 
the particle oscillates around one of the two new minima. Clearly in this case we find the perfect 
agreement between the exact and the TWA solutions for considerably long times. 

For a sudden quench 5 — > oo the classical equations of motion governing the dynamics deter- 
mined by the Hamiltonian (|117p can be solved explicitly: 

x(t) = V2^cs (VT+\(u>t + 13)1 + T _^_) , (119) 

where cs(u|m) is the elliptic function and a and (5 are the free parameters determined by the initial 
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Time 



FIG. 2 Same as in Fig. [TJbut for a sudden change of the parameter k from 1 to —1. It is clear that the 
TWA accurately describes the first oscillation of (x 2 (t)) but then starts to deviate from the exact result. 

conditions. If the parameter A in Eq. (|117p is small then at short times the quartic term in the 
Hamiltonian is irrelevant and Eq. (|119p simplifies to 

x(t) ps y/2a sinh(/3 + ut) (120) 

so that 

xo ps v / 2asinh/3, po ps \^2auj/m cosh(/3). (121) 

These relations allow one to express the unknown constants a and (3 through initial conditions and 
then average an observable of interest over these initial conditions with the Wigner function. In 
particular, at short times when the approximation f| 1 20 [) still holds (or when x 2 <C 1/A) we find 
that 

(x 2 ) « — — cosh2o;t (122) 

From this expression we can estimate the time it takes the particle to reach one of the minima of 
the mexican hat potential 



2u 



2mw 



(123) 
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FIG. 3 Same as in Fig. [2] The additional blue line represents the first quantum correction evaluated 
according to Eq. (|53|) with At = 0.25. Despite this value of the time step is quite large, we checked that 
the results are insensitive to this particular choice. It is clear that the quantum correction extends the 
agreement of TWA with the exact solution to a longer time. 

This result can be of course obtained by more elementary methods. Using properties of the elliptic 
functions one can analyze corrections to Eq. (|122p . Our goal here is, however, to understand the 
limits of applicability of the semiclassical approach and analyze the lowest quantum corrections to 
the dynamics. For this purpose we will proceed with numerical results. 

In Fig. [2] we plot the results of the semiclassical (TWA) and exact calculations for the variance 
of the displacement of the particle after a sudden quench. It is clear that TWA describes very 
well the initial spread of the wave packet representing the particle and its following collapse. 
However, at later times the exact and semiclassical curves start to depart from each other and 
the quantitative agreement is lost. In Fig. [3] we show the same plot on a shorter time scale with 
an additional line representing the first quantum correction evaluated according to Eq. (|53p . This 
correction extends the validity of TWA, however, always remaining relatively small and thus failing 
to indicate where TWA breaks down. As we will see in the next example the situation can be quite 
different and the quantum correction can accurately predict where TWA stops being valid. Note 
that for evaluating the correction formally one needs to take the limit At — > 0. However, in 
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practice we used a relatively big time step At = 0.25 and yet got accurate results, which do not 
change if At becomes smaller. From this simple example it is evident that the quantum diffusion 
gives very robust results similarly to the classical Langevin type diffusion and thus can be used for 
practical calculations. We note that the formal evaluation of the response according to Eq. (|52p 
gives very poorly converging integrals because of the unstable classical dynamics near x = and 
thus it is not suitable for numerical purposes. 



C. Sine-Gordon model. 



Further increasing the complexity of the system in this section we will show how the formalism 
can be applied t o the sine-Gordon (SG) model, which is described by the following Hamilto- 



nian ( Giamarchi 



200J): 



n 



dx 
T 



h 2 {x) + [d x <t>(x)Y - 2V cos(/30(x)) 



(124) 



where h(x) and <ft(x) are conjugate variables with canonical commutation relations (like coordinate 
and momentum) and is the parameter, which characterizes the strength of quantum fluctuations 
and plays the role of H. This model has numerous applications in different areas of p hysics. 



19971 ). The 



Both classical and quantum SG models are integrable (jLukvanov and Zamolodchikovl . 
discrete version of the SG model, where dp(j> — > ( frj+i — 4>j, is known under the name of the Frenkel- 



Kontorova model rtChaikin and Lubenskv 



1 



T + 2a 2 



19951 ). Its Hamiltonian is 

2 -Vcos(0^] 



'Jj+ij 



(125) 



Here a is the lattice constant with a — > corresponding to the continuous SG model. The model 
above has identical low-energy properties as the SG model. We will use the Frenkel-Kontorova 
model in actual numerical simulations and will set a = 1. 

If the SG model represents commensurate boson s in an optica l lattic e potential then = 



I^J'K /K, where K is the Luttinger-liquid parameter ( Biichler et all . I2003) 3 . For < 2V 27T the 



cosine potential is relevant in equilibrium and for > 2\/2ir it is irrelevant (jGiamarchi 



200J). This 



means that in the former case even for t 



spectrum (IDonohue and Giamarchi 



l e infinitesimal value of V the system opens a gap in the 



2001 



Mathey et al. 



200a). In this regime the structure of the 



3 There is some controversy in the literature in defining either the L uttinser-Liquid parameter or its invers e enters 
the expression for j3. We use the convention discussed in e.g. Refs. dCazalillal . l2004j ; |Pe Grandi et a/.ll2008l ). where 
large K corresponds to weak interactions 
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spectrum of the SG model is quite rich. Thus at (3 < 2yfir the elementary excitations, solitons, can 
bind together to form new qua siparticles B\ breathe rs. Then at (3 < \[2~n B\ breathers can bind to 



form £>2 breathers and so on (jGritsev et al 



2007a 



bj). Details of the spectrum of the excitations 



are important for understanding non-equilibrium dynamics. In particular, it was shown that for 
(3 < quench in the SG model can lead to damped or u ndamped oscillations at frequencies 



corresponding to the breather's masses (jGritsev et al. 



2007al ). For large /3, where the cosine term 



is irrelevant and the spectrum is gapless, one can expect strong damping of the oscillations or no 
oscillatory behavior at all. 

For the purpose of this paper it is important to realize that f3 plays the role of the Planck's 
constant. This can be seeing, for example, from the rescaling (/>//? and noticing that 1/(3 
becomes the saddle point parameter of the quantum evolution operator. Thus for f3 <C 1 we expect 
that the semiclassical description is valid and for f3 > 1 it fails. To simulate the time evolution in the 
SG model we chose a particular dependence of the coupling V on time: V(t) = 0.1 tanh(0.2£). For a 



syste m of commensurate bosons this corresponds to turning on an optical lattice (|De Grandi et al 



2008). We will be interested in the expectation value of cos(/3(/>) as a function of time. On general 



grounds we anticipate a large response for small f3, where the potential is strongly relevant, and 
conversely a small response at large (3. 

First let us assume that the system is initially prepared in the ground state. The latter corre- 
sponds to the product of ground states of harmonic oscillators corresponding to different momenta 
q. Thus the Wigner function of this state is a product of Gaussians. In particular, for the Frenkel- 
Kontorova model with a = 1 we have 



W(<P° q ,n° q ) = l[exp 



,012 

"' 2a 2 Jn° n \ 2 



2a\ 



(126) 



where o q = 1/ y / 2u q , ui q = 2sin(g , /2), and q = 2i:n/L with n = 0, 1, ..L — 1 (L is the system size). 
We note that the zero momentum mode is characterized by n q 
distribution we can generate the initial values for (fp- and in real space: 



and random <p q . From this 



(127) 



q q 

Within TWA we need to solve the classical equations motion: 



dt 2 



+ i _ 1 -2^-)- /3F(t)sin(/%) 



(128) 



n° distributed ac- 



subject to the random initial conditions <pj(t = 0) = 0°, 4>j{t = 0) 
cording to Eqs. (|126p and (|127|) . The first quantum correction to the semiclassical approx- 
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imation can be evaluated according to Eq. (|53p . where the third derivative of the potential 

Before showing numerical results we would like to emphasize the difference between expanding 
evolution in quantum fluctuations (j3) and in the strength of interaction V. Note that since for 
noninter acting problems TWA is exact the expansion in j3 is also the expansion in V . However 
these two expansions are very different. The expansion in f3 requires solving nonlinear classical 
equations of motion taking into account V at the classical level in all orders of perturbation theory. 
Conversely linear response (and higher order corrections in V) assumes that the potential is a weak 
perturbation but treats all orders in /3 exactly. It is clear that the perturbative expansion in V 
should work better at large values of f3 where the cosine potential is less relevant. Conversely the 
expansion in (3 should work better when f3 is small and the cosine potential is strongly relevant . 



The expectation value of cos(/3(/>) within the linear response theory is given by (jAbrikosov et al 



19651 ): 



(cos /90(t)) ~ i / dt dxV(t')( cos(/30(x,O),cos(/30(O,f)) 



(129) 



where the index "0" indicates that the averaging is taken over the initial state with the evolution 
given by the unperturbed Hamiltonian (without the cosine term). The correlation functio n above 



is related to the imaginary part of the retarded correlation function (jAbrikosov et al 



19651 ). which 



can in tu rn be evaluated b y analytic continuation of the corresponding imaginary time correlation 



function (jGiam archil . 



2004) to the real time axis. Then we find: 



(cos « / dry(t-r)^exp 



/3 2 ^-4 1 — cos(gj) cos(cl> 9 t) 



p \ ^ 
2L ^ 



q=0 



sm 



/3 2 ^4 cos(gj) sin(a; (? r) 



2L ^ 



g=0 



UJ (I 



(130) 



We note that the q = mode here has to be included in the expression above with the 0/0 
uncertainty evaluated according to the L 'Hospital's rule. 

We can also evaluate the analytical expression for the linear response within the semiclassical 
truncated Wigner approach: 



(cos/3cHt)) 4 



drV (t — t) ^2 ex P 



2L ^ 



COs(qj) COs(u! q T) 



q=0 



2 



L-1 



P \ ^ 

2L ^ 



cos(qj) sin(ui q r) 



g=0 



(131) 



The expressions (|130p and (|13ip almost coincide except for the last multiplier. Clearly in the limit of 
small (3 the argument of the sine function in the exact expression is small and the difference between 



43 



the two results vanishes. Therefore within the linear response the TWA result is accurate up to 
the terms 0(/3 6 ). Comparison between Eqs. (I130p and (I13ip also shows that TWA gives accurate 
results at short times even if f3 is not small. This is in accord with our general expectations that 
TWA is asymptotically exact at small t. Whether TWA always remains accurate or diverges at 
long times is on the other hand highly non-universal. For example if the sum over momenta q in 
the argument of the sine function in Eq. (|130p is bounded at all times, which would be the case 
in higher dimensions, then the difference between Eqs. (|130p and (|13ip remains small at all times 
provided that /3 <C 1. One can expect that in this situation TWA remains accurate for infinitely 
long times even beyond the linear response approximation. 

Next we compare the results for the expectation value of cos(/30) within TWA and the additional 
quantum correction with the full linear response calculations for different values of /3. As we 
mentioned above we take V(t) = 0.1 tanh(0.2i) so that the potential always remains small and the 
linear response is justified at least for relatively short times. Thus we can use the linear response 
analysis to determine the accuracy of the semiclassical approach and the quantum corrections at 
short times. 



0.10 




FIG. 4 Time dependence of (cos(/30)) for the discretized SG model described by the Hamiltonian (|1 24[) with 
the system size L = 20, f3 — 2y/ir, and the cosine potential turning on as V(t) = 0.1 tanh(0.2). The three 
plots correspond to the TWA, TWA with the first quantum correction, and the linear response. 
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In Fig. H] we show time dependence of (cos(/3(/>)) for /3 = and L = 20. The three plots 
correspond to the semiclassical approximation, the first quantum correction evaluated according 
to Eq. (153p with At = 0.2 (see Appendix [A] for additional details), and the linear response. While 
the cosine potential is relevant in this regime the linear response is remains valid for sufficiently 
long times and thus gives an accurate prediction for the exact result. It is clear from the plot that 
TWA is valid only at very short times and then rapidly breaks down. The first quantum correction 
significantly extends the TWA approach and unambiguously shows where TWA breaks down. 




FIG. 5 Same as in Fig. (0} but for (5 = V2n~. 

The situation changes dramatically if we consider smaller values of f3, i.e. if we go closer to the 
classical limit. Thus in Fig. [5] we plot (cos (5(f)) for (5 = V2tt, i.e. a factor of y/2 smaller than in 
Fig. |U Despite this is a relatively small change in the parameter f3 the results are quite different. 
First the linear response is clearly valid for much shorter times t < 5. After that the linear response 
calculation is not reliable leading to a spurious divergence of (cos /3</>) . Second the accuracy and the 
time domain of validity of the semiclassical approximation is significantly improved compared to 
the previous case. The first quantum correction again clearly extends the domain of applicability 
of the semiclassical approach. And finally in Fig. [6] we plot the same dependences for even smaller 
/3 = y/n. Now the linear response result shows a very good agreement with the TWA calculation for 
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FIG. 6 Same as in Fig. (gj} but for /3 



relatively short times t < 5. The first quantum correction makes this agreement almost ideal. At 
longer times the linear response clearly breaks down, however the first quantum correction to the 
semiclassical result continues to be small. Thus we can expect that the semiclassical approximation 
is very accurate here. It is interesting that by changing the parameter (3 (which plays the role of 
h) by only a factor of two we expand the regime of validity of TWA by at least two orders of 
magnitude. 

Instead of starting from the ground state one can consider a situation, where the system is 
initially prepared at finite temperature T. The Wigner function for the thermal density matrix 
can be straightforwardly evaluated using Eq. (|148|) (see also the supplementary information of 



Ref. rtPolkovnikov and Gritsev 



2003)): 



W(4> q ,n° q ) 



JJexp 

q 



m 2 



where r q = coth(a; g / (2T)). In the limit T -C oj q we have r q 
result (]126p . In the opposite limit T S> 0J q we get 



(132) 

1 and we recover the zero temperature 



W(</>° q ,n° q ) 



JJexp 

i 



^4 + Kl 2 

2T 



(133) 
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which is nothing but the Boltzmann's distribution. We note that for any nonzero temperature the 
distribution of the zero momentum mode (q = 0) is always Boltzmann's because there is no zero 
point energy associated with this mode, which simply describes translation of the whole chain. 
It is interesting to observe that the uncertainty between conjugate variables 4>° q and n° q at finite 
temperatures satisfies: 

Thus we recover the minimum uncertainty relation at T — > and <5|<^|J|<5|ng| ~ 2T 2 /K| 2 > 1 at 
large T. 

It is also straightforward to generalize Eq. (|130p to the case of nonzero temperatures: 



(cos/3<£(t)> « f drV{t-T) 
Jo 




.(135) 



FIG. 7 Same as in Fig. (j4j but for /3 = V2n and initial temperature T — 0.2. 

The results of the TWA simulations, the first quantum correction as well as the linear response 
calculations for T = 0.2 and /3 = \/2tt are plotted in Fig. [71 Perhaps a somewhat unexpected 
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feature of the simulations is that even at finite initial temperatures the semiclassical approximation 
still shows significant deviations from the exact result. As in the zero temperature case the first 
quantum correction considerably extends the time domain of the applicability of the semiclassical 
approximation. We want to stress again that the small parameter of the theory is j3. In simulations 
we used f3 = \p2/n corresponding to the regime of relatively strong quantum fluctuations. We 
intentionally used this regime to emphasize the significance of the quantum correction. For smaller 
value of j3 the accuracy of TWA improves dramatically like in Fig. [6j 



D. Dynamics in the coherent state representation. 

In this section we will address various problems in the coherent state phase space. We will review 
both examples discussed in some earlier works and present a couple of new ones. Most of the exam- 
ples here have direct relevance to the cold atom experiments. We want to emphasize that the main 
goal of this section is to verify accuracy of the formalism discussed here. So some of the problems 



tally relevant situatio ns in cold atom systems can be found in the review 



in more recent works ( 


Fane: et al. 


, 2009 


; HiDolito and Polkovnikov 


2009; 


Mathev and Polkovnikov 


, 2 


009 


; £ 


>au et al.. 


2009 


)■ 



20101 : 



1 

(Blakic et al, 


2008 


) and 


Martin and Ruostekoski. 



1. Collapse of a coherent state 



The first example we will briefly review here deals with the evolution of the initial coherent state 
with on averag e N particles per si te subject to the Hamiltonian containing only the interaction 
term (see Ref. (|Polkovnikovi . l2003al ) for more details): 

U 



-^tp(^tp - 1). 



(136) 



This problem is closely related the collapse-revival experiment by M. Greiner et. al. (jGreiner et al 



2002). Because the problem does not have a kinetic term it can be easil y solved ana. 



particular, the expectation value of the annihilation operator is given by (jPolkovnikov 



ytically . In 



2003al ): 



($(t)) = v^exp [N(e 



-iUt 



1) 



(137) 



This result can be obtained by expanding the coherent state in the Fock basis: 



-N/2 



n) 



(138) 
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and then trivially propagating this state in time. The classical limit is achieved here taking the 
limit N — > oo, U — > while keeping the product UN = A fixed. The exact evolution shows a rapid 




FIG. 8 Dependence of the real part of (ip) on time for the initial coherent state evolving according to the 
Hamiltonian (|136|) . The interaction is taken to be U = 0.5 and the average number of particles in the initial 
state is N = 4. 



collapse of the coherent state at time t c ~ 1/UN = 1/A and then the consequent revival at time 
t r ~ 1/U = N/X (see Fig. [8]). In the classical limit t r — > oo thus the revival is a purely quantum 
phenomenon. 

Next we solve the problem using TWA and quantum corrections. For doing this we need the 
Weyl symbol for the Hamiltonian (|136p : 



UwW, , *l>) 



which yields the following classical equation of motion 

•W) rrn,UM2 



at 



u(m)\-mt)- 



(139) 



(140) 



This equation can be trivially solved using that \ip(t)\ 2 = \ipo\ 2 is the integral of motion. The 



solution s 



rould be supplemented by random initia 



function (jGardiner and Zollerl . 



2004 



Polkovnikovl . 



condi tions distributed according to the Wigner 



2003aJ ): 

W{Ml) = 2exp[-2|Vo - ^| 2 ] 



(141) 
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Using the explicit analytic solution of Eq. (|140p and the Wigner function above one can ana- 
lytically cal_cj : dj^e_l2pth_the TWA result (ifto(t)) and the first (and higher) quantum correction 



{?pi(t)) (IPorkovnikovi . l2003ah : 

ip (t) = y/N exp 



iUNt 
' 1 + iUt/2 



exp[iUt] 



U 2 t 2 

^(t)=^(t)—\l 



iUNt 



(l + iUt/2) 2 ' 



iUt 



(142) 



(143) 



4 V~ 3(l + iUt/2) 2 3(l + iUt/2) t 

It can be explicitly checked that the Taylor expansion of the exact result for {tp(t)) in powers of 
1/N (at a fixed A = UN) coincides with the similar expansion of the TWA result ipo(t) up to 
the terms of the order of 1/N 2 and with the expansion of ipo(t) + ipi(t) up to the terms of the 
order of 1/N 4 . To illustrate these dependences we plot the real part of the expectation value of 
(tp(t)) versus time in Fig. [8] for the parameters U = 0.5 and N = 4 (corresponding to A = 2). 
It is clear that at short times TWA gives a very good approximation to the exact result and the 
next quantum correction improves it even further. However, both TWA and the correction fail to 
describe the revival phenomena which have a purely quantum origin. From Eq. (|137p one can see 
that in order to correctly describe the revival one needs to include exponentially large number of 
terms. Mathematically this happens because the functions of the type exp[exp[ix]] have very bad 
convergence of the Taylor expansion at large x. 



2. Evolution towards the Schrodinger cat state. 



Next we will consider an example largely following Ref. (|Polkovnikovi . l2003bl ) dealing with 
dynamics of the system of bosons in a double well potential. We assume that the only two lowest 
single-particle (vibratio nal) states are relevan t so the system can be described by the Bose-Hubbard 
Hamiltonian (see Ref. (ITrotzky et all . l2008l ) for discussion of applicability of the Bose-Hubbard 
model to cold atoms in optical lattices): 



n 



(144) 



where j = 1,2, J is the tunneling matrix element and U is the interaction constant. We assume 
that initially the system of 2iV particles is prepared in the noninteracting ground state and then 
slowly attractive interactions are ramped up. We note that alternatively on e can initially prepar e 



the system in the most excited state and then turn on repulsive interactions (|Polkovnikovi . 



2003b!) . 



Beyond some critical strength the classical symmetric equilibrium, where both wells are equally 
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populated, becomes unstable and the system evolves towards the "Schrodinger cat" state. As in 
the previous example the classical limit is achieved when the average number of particles per site N 
becomes large while the product UN is kept fixed. The corresponding classical (Gross-Pitaevskii) 
equations read 

.dMt) 



dt 
.dMt) 

'' dt 



-j^ 2 (t) + u(t)(\Mt)\ 2 -i) Mt) 
-jMt) + u{t){\Mt)\ 2 -i) Mt)- 



(145) 
(146) 



As in Eq. (j!40j) extra "—1" in Eqs. (I145p and (jl46H appear because of the Weyl ordering of the 
Hamiltonian. Since these terms only result in the overall phase of the wave function they can be 
dropped. We assume that initially the noninteracting system of 2N particles is prepared in the 
ground state and then one slowly introduces an attractive interaction: 

1 tanh(cft) 



U(t) 



(147) 



N I -St ' 

where 5 is the parameter controlling the ramp speed. We have chosen such dependence of U{t) 
so that at time t = 1/5 the ratio of the interaction to the tunneling becomes infinite. Physi- 
cally this corresponds to linearly shutting down the tunneling between the two sites since only 
the ratio J/U is physically important. The parameter 5 in Eq. (j!47H controls the speed of the 
ramp. We note that this problem is very closely related to the motion of the particle in the Mex- 
ican h at potential if we use number imbalance n = t/>|t/>i — i^\i^2 as an effect i ve co ordinate (see 



Refs. (jHipolito and Polkovnikovl . l2010l : 



Polkovnikov 



2003a 



Polkovnikov et al. 



2002 ) for details). 



At sufficiently weak interactions the particles evenly distributed between the two wells and the ef- 
fective potential has a minimum at n = 0. As interactions exceed some critical value U c ~ J/N the 
equilibrium n = becomes unstable and the particles prefer to concentrate in one of the wells. For 
large number of particles the total number constraint is not important and the initial s tate is almost 



equiva lent to the product of coherent states, for which the Wigner function reads (|Polkovnikovl . 



2003TJ ): 



W(ipj,ip1j) = I}2exp [-2|Vv - y/Nexp(i 



(148) 



The overall site independent phase here is arbitrary and we choose = 0. For the initial state 
where t he number o f boson s is fixed the Wigner function is also readily obtained (see Eq. (43) 



in Ref. ( Polkovnikov 



2003bl )). however it contains highly oscillating component and is thus much 



harder to deal with. In large systems the difference between coherent state and fixed number 
symmetric state is not important. In Fig. [9] we show the dependence of the number variance on 
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FIG. 9 Dependence of the normalized number variance Sn = \J (h 2 ) /2N on the interaction strength. The 
number of particles per cite is N = 32 and the interaction depends on time according to Eq. (| 147[) . The 
solid black line represents the exact solution and t he dashed blue line is the result of the semiclassical 
approximation (TWA). This plot is taken from Ref. (|Polkovnikovl l2003bf ) . 



the interaction strength for three different ramp rates 5. It is clear that the semiclassical approach 
gives very accurate predictions for slowest and fastest rates with small deviations occurring only 
at the intermediate rate. 



3. Turning on interactions in a system of lattice bosons. 

Let us now consider a process of heating bosonic atoms due to non-adiabatic changes of the 
coupling in time. As in Sec. IIV.D.2I we assume that the system is described by the Bose-Hubbard 
model but with more than two sites: j = 1, . . . , M. For simplicity we chose periodic boundary 
conditions. We analyze a process where one starts from the ground state in the noninteracting 
regime and then ramps on now repulsive interactions according to 

U(t) = U tanh(St). (149) 

The classical Gross-Pitaevskii equations of motion are straightforward multi-site generalization 
of Eqs. ()146p . To avoid dealing with highly oscillating Wigner function, describing the initial 
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conditions, it is convenient to start from the coherent state rather than the state with fixed number 
of particles. Then the Wigner function is given by Eq. (I148p . To see the effects of quantum 
corrections we intentionally choose parameters correspondi ng to strong quantum fluctuatio ns: J = 



20081 ) we will 



1, C/o = 1 and on average one particle per site. As in Ref. (jPolkovnikov and Gritsevl . 
analyze the energy density (or the energy per site) in the system. 

First we consider a small periodic one-dimensional chain consisting from eight sites so that we 
can do comparisons with exact results obtained by directly solving the Schrodinger equation. In 
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FIG. 10 Dependence of the energy density on the product St (see Eq. (|149jl ) for a periodic Bose-Hubbard 
chain consisting of eight sites with on average one particle per site. The hopping amplitude J = 1 and 
interaction changes in time according to Eq. (| 149|) with Uq = 1. Initially the system is assumed to be in the 
coherent zero-momentum state. The two sets of curves correspond to two different ramping rates <5 = 0.5 
(solid symbols and a solid line) and S = 2.5 (open symbols and a dashed line). Red circles correspond to 
TWA and blue squares include the first quantum correction evaluated according to Eq. ([75]). Black lines 
represent exact results obtained by direct solution of the Schrodinger equation. 



Fig. [10] we show comparison between exact results as well as TWA and TWA + first correction 
(evaluated according to Eq. (|75p ) for the two values of ramping rate: 5 = 0.5 and 8 = 2.5. We 
discuss details of numerical implementation of the quantum correction in the Appendix El It is 
evident that TWA gives rather accurate results even though the system ends up in the strongly 
interacting regime. The first quantum correction significantly improves the accuracy of TWA 
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making the results practically indistinguishable from the exact ones. This example highlights 
again that the quantum correction can be used not only to significantly improve the accuracy of 
TWA but also as a method to find where TWA starts to deviate from the exact result. 

Next we will consider the same setup in a two-dimensional system of 32 x 32 sites. This system 
is far beyond the limits of exact quantum simulations. Moreover this parameter regime corresponds 
to relatively strong interactions and significant heating However, it can be readily simulated using 




FIG. 11 Same as in Fig. [T0]but for a two-dimensional 32 x 32 lattice. We use 5 = 1, J = 1, on average one 
particle per site. The solid lines correspond to Uq — 2 while the dashed lines do to Uq = 1. As before red 
circles represent TWA and blue squares describe TWA + first quantum correction. 

the phase space methods discussed here. In Fig. [11] we plot energy per site for the ramp with 5 = 1 
and two values of the final interaction Uq = 2 (solid lines) and Uq = 1 (dashed lines). We also use 
J = 1 and the unit filling. As expected TWA is more accurate for smaller interactions, which is 
reflected in a smaller difference between the blue and the red graphs. 

E. Thermalization within the Bose-Hubbard model. 

Let us now analyze the applicability of TWA to a very important issue of thermalization. We 
will again focus on the Bose-Hubbard model discussed above. As a specific situation we assume 
that initially all atoms are placed into one of the sites and then they are allowed to spread over. 
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Since the system is nonintegrable one can expect that it wi 



there is a famous example of the Fermi-Pasta-Ulam problem (jFermi et al. 



1 thermalize at long tim es. However 



May 19551 ). which shows 



that the thermalization is not necessarily guaranteed even on nonintegrable systems. Overall the 



issue of thermalization in c 



tention (Kollath et al. 



2007; 



osed Hami l tonia n systems at t racted recently a lot of theoretical at- 



Reimannl . 



2008; 



Rigoj, 



2009; 



Rigol et al. 



20081 ). Because of lack of 



available tools to address dynamics in nonintegrable systems one either has to rely o n exact simula- 



tions like in Ref. (jRigol et al. 



2008) or on DMRG and related methods like in Ref. (jKollath et al. 
20071 ). The former can deal only with small system sizes, while the latter have long time limita- 
tions and generally applicable only to systems close to equilibrium. From this prospect phase space 
methods, in particular TWA, provide a very natural framework to study thermalization dynamics. 
In particular, TWA and further corrections give asymptotically accurate results at short times. At 
long times, if the system thermalizes, quantum fluctuations become subdominant if the temperature 
characterizing the steady state is not small. Then quantum fluctuations should remain subdom- 
inant and TWA results are expected to continue to be valid. While this discussion is somewhat 
speculative and requires further thorough analysis we will show here that these simple ideas are valid 



for our s imple toy example. Therma. 



i zation within the phase space approach was also analyzed 



in Re 


fs. ( 


Bistritzer and Altman. 


2007: 


HiDolito and Polkovnikov. 


2010; 


Mathev and Polkovnikov. 


2009; 


Polkovnikov and Gritsev 




2008 


i 


Polkovnikov and Wang. 


2004: 


Wright et al. 


2008 


) and oth- 



ers. We assume that the system is initially prepared in the coherent state in a single site with 
an average number of atoms N. In Fig. [12] we show the dependence of the occupancy of the zero 
momentum state for such initial conditions for a periodic chain of six sites. 

The classical simulations show very little sign of thermalization. In this sense this problem is 
very analogous to the Fermi-Pasta-Ulam problem with the main difference that the initial state is 
localized in the coordinate rather than momentum state. Once we include quantum fluctuations 
already at the level of TWA we clearly see that the system approaches a steady state. Comparing 
these results with those of the exact diagonalization we see that TWA gives extremely accurate 
predictions. Obviously for the situation shown the quantum corrections are negligible (and this 
indeed can be checked numerically). It can be checked that at longer times the exact quantum 
simulations show revival of the initial state due to a small system size. As in earlier examples TWA 
misses that quantum revival. We expect that these revivals disappear as we go to larger system 
sizes. Since the goal of this paper is to review quantum dynamics using the phase space methods, 
we will set aside the question whether the steady state seen in Fig. [12] corresponds to the thermal 
equilibrium and address it in a separate publication. 
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FIG. 12 Dependence of the occupancy of the zero momentum state on time for the initial state where all 
particles populate a coherent state in a single site. The total number of sites is 6, the average number 
of particles is N = 18, the interaction strength is U = 0.15, and the tunneling is J = 1. The blue curve 
describes the classical (Gross-Pitaevskii) simulations, which show no sign of thermalization. The black curve 
is the result of exact simulations and the red curve corresponds to TWA. 

We can easily extend TWA and classical GP-simulations to larger arrays, where the exact 
integration of the Schrodinger equation is no longer possible. In Fig. [13] we show the evolution 
of the zero momentum occupancy for the initial array of coherent states having the same phase 
and occupying each tenth site. (We start from coherent states because they correspond to the well 
defined classical initial conditions). TWA again clearly indicates very fast thermalization while the 
classical curve shows non ergodic behavior. 

F. Spin dynamics in a linearly changing magnetic field: multi-level Landau-Zener problem. 

In this section we will illustrate how phase space methods can be applied to another exactly 
solvable problem. Namely we will consider a motion of an arbitrary spin S in a linearly changing 
magnetic field: 



H = 2h z {t)s z + 2gs x , 



(150) 
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FIG. 13 Same as in Fig. [12] but for an array of 60 sites. The initial state is such that each tenth site is in 
a coherent state with 18 particles on average and the same phase. The two curves show the results of the 
semiclassical (TWA) and the classical GP simulations. The TWA apparently predicts rapid thermalization, 
while the classical evolution is clearly non-ergodic. 

where h z (t) = St. We assume that the system is initially prepared in some state or superposition 
of states at t = —to and we will be interested in finding expectation values of various observables 
at t = to, where to is large so that h z (to) S> g. 

As we discussed in Sec. IIII.CI one can map the time evolution of noninter acting spins to the 
evolution of noninteracting bosons using the Schwinger representation. Therefore TWA is exact in 
this case. Using Eqs. (|84h the Hamiltonian (I150p becomes: 

M = h z {t)(a ] a - /3 f /3) + g(a ] '/3 + /3 f a). (151) 

The classical equations of motion corresponding to this Hamiltonian are 

(Lol 

i— = 5t a + g(3, (152) 
at 

dB 

iJL =ga - StB. (153) 
at 

These equations should be supplemented by the initial conditions distributed according to the 
Wigner transform of the initial wave function (or the density matrix) . 
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Note that Eqs. (|152p and (|153p map exactly to the equations describing the conventional Landau- 
Zener problem. Then the evolution can be described by a unitary 2x2 matrix: 



Ta + R&Pa, /?oc = -JfcT^ao + T/3 



where (see Refs. (jKeeling and Gurarid . 120081 : IZenerl . Il932l )) 



T = e"^ 7 , R = yl — T 2 , = 7 [ln( 7 ) - 1] - 2-fln(V26T), 



(154) 



(155) 



and 7 = g 2 /(25) is the Landau-Zener parameter. 

Finally we need to reexpress different spin components at t = 00 through the initial values of 
the spins. Thus we have 



4, = (T 2 - fl ^foo^A) + Q5Aji£re ^ + cofiKre-i* 

= (T 2 - i? 2 )sg + 2i?T cos(0)sg - 2i?T sin(0)4, (156) 
= -2RT cos{$)s z Q + (T 2 - i? 2 cos(2</>))s§ + # 2 sin(20)sg, 
= 2i?Tsin(^.)sg + i? 2 sm(2<£)sg + (T 2 + # 2 cos(20))sg. 

Now using these expressions and the results of Sec. IHI.CI we can compute expectation values of 
various quantities. For simplicity we assume that initially the spin is prepared in the eigen state 
of s z . On the language of bosons this corresponds to the Fock state \S — n,n), where a particular 
value of n corresponds to the state with s z = S — n. Let us give a few specific examples: 

<rj = (T 2 - R 2 )s z 

((s^) 2 ) = [T 4 + R 4 - AT 2 R 2 ] (s z ) 2 + 2T 2 R 2 s{s + 1), 

+ Zo%o) = ^RT(T 2 - R 2 ) cos(0) [s(s + 1) - 3<(s§) 2 >] . (157) 

Note that for the conventional Landau-Zener problem corresponding to the spin s = 1/2 the last 
two equations become trivial: ((s^,) 2 ) = 1 and (s^s^ + s^s^) = 0. But for larger values of 
the spin these correlation functions are nontrivial with e.g. (s^s^ + s^s^) being an oscillating 
function of the total time to an d the Landau-Zener parameter 7. 



G. Dicke model. 



As a final illustration of the phase space approach to quantum dynamics we will consider the 
Dicke model (or more accurately the Jaynes-Cummings model), which maps to the bosonic mode 
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interacting with a large spin: 



n 



-\&b + -jL= ( Ps~ + hfr 



(158) 



where b is the bosonic annihila tion operator and 5 + , 5 are the raising and lowering spin operators 



(see Ref . (jAltland et all 120091 ) for additional discussions) . Here 5 is the value of the spin and the 
classical limit is achieved at 5 — > oo. We will specifically consider a setup where one starts at A 
large and negative in the state where bosonic mode is unoccupied and the spin is polarized along the 



z-axis. This problem was analyzed in Ref. (jAltland et al. 



2OO9I ) using various methods including 



TWA. Here we will repeat the details of TWA analysis to this problem because it illustrates the 
power of phase space methods to describe all aspects of quantum dynamics. This problem also 
illustrates the possibility of mixing different phase space representations together, namely the 
coherent state phase space describing the bosonic mode and the 50(3) vector describing the spin 
degree of freedom. This p roblem is particularly important also in view that using the Anderson 

,581) describes the system of degenerate fermions 
2OO9I ) . Thus the dynamics in this case represents 



mapping ([Anderson . 



19581 ) the Hamiltonian 



interacting with a bosonic mode (jAltland et al. 
a si mplified versi on of the dynamics through the BEC-BCS crossover in cold atoms (see E.g. 



Ref. ((RegaJ, 



20061 )). Under this mapping in the classical limit (which is equivalent to the mean- 



field approximation in the BCS model) the dynamics of t he system can b e ful 



syste m of classical equations of motion analyzed in Refs. (jAndreev et al 



2004 



y mapped to the 



Barankov et al. 



2004) for a more general case with the dispersion. TWA and corrections show the way to go beyond 



the purely classical limit and indicate that the whole approach of this paper can be extended to 
interacting fermionic systems at least in certain situations. We note that the model (1158j) also 
appears in various other contexts in optics, condensed matter and atomic physics. 

We choose a specific process where A changes linearly in time: A = 25t. Initially at t — > —00 
there are no bosons and the spin points along the z direction: (S z ) = 5. As the time increases 
the bosonic mode becomes resonant with the spin forcing it to rotate. Eventually when time 
becomes large and positive the bosonic mode becomes off resonant again. For 5 = 1/2 this 
problem is identical to the Landau-Zener problem. For larger values of the spin it becomes a 
particular generalization of the Landau-Zener problem to multiple energy levels. However, unlike 
in the previous example, this generalizat ion results in an interacting quantum prob l em w hich has 



very nontrivial adiabatic (5 —> 0) limit (jAltland and Gurarie . 



2008 



Altland et al. 



20091 ). Using 



Eqs. (|95p and (|148p with N = we find that the Wigner distribution of the initial conditions for 



59 



large S can be well approximated by 

W(b*,b,S x ,S y ,S z ) « 2e- 2 W 2 -^e-( s * +s 'tV s 5(S z - S). (159) 
The classical equations of motion corresponding to the Hamiltonian (jl58R read: 

* = -^b x s, (iao) 

where B = ($lb, 0) is the effective magnetic field acting on the spin. 




FIG. 14 Average relative number of bosons rt& = (wb)/2S created during the process as a function of the 
rate of change of the chemical potential 5. Initially the bosonic level is empty and the spin points in the 
positive z direction. Shown are the semiclassical and exact results for two values of spin S = 32 and S = 64. 
The inset shows the difference between exact and semiclassical results m ultiplied by a f actor of 100. The 



parameters of the model for the simulations are identical to those in Ref. (jAltland et al 



2009) 



In the purely classical limit the initial conditions are b = 0, S z = S, S x = S y = 0. It is easy 



to see that the solution b(t) = 0, S z (t) = S satisfies both the classical equations of motion and 
the initial conditio ns. It is thus necess ary to include quantum fluctuations to see the nontrivial 



dynamics. In Ref. (jAltland et al 



20091 ) it was shown that using conservation laws Eqs. (|160p can 
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be simplified further and be reduced to a single second order differential equati on. The latter 



in tu rn can be approximately solved using either methods of adiabatic in variants (lAltland 



et al 



200 



2009) or more accurately using the exact solution of the Painleve equation (jltin and Torma . 
This solution indicates a very nontrivial adiabatic limit for large values S. Since our goal here is 
to understand the validity of the TWA we will stop short of discussing details of this solution and 
focus on the comparison between the results of exact simulations and TWA. In Fig. [14] we show 
the average relative number of created bosons n b = (¥b)/2S during the process described above 
as a function of the rate 5. Slow rates correspond to the adiabatic limit and thus to nearly 100% 
conversion efficiency (nearly maximum possible number of bosons in the final state). We plot the 
exact and TWA results for two values of spin S = 32 and S = 64. The accuracy of the semiclassical 
calculation in the whole range of 8 does not exceed 1% for S = 32 and is about a factor of two less 
for S = 64 (see the inset in Fig. [T4l) . 




FIG. 15 Fluctuations of the relative number of bosons Snt, = y ((b^b) 2 ) — (wb) 2 /2S as a fun ction of 5. The 



param eters are the same as in Fig. H! with total spin S = 32. This plot is taken from Ref. (jAltland et al. 
20091) . 



We can also check that the semiclassical approximation is accurate for other observables. For 
example in Fig. [15]we plot fluctuations of the relative number of bosons at t — > oo as a function of 5. 
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The number fluctuations are small both in the limit of large 8 where almost no bosons are created 
and in the limit of small 8, where there is almost 100% conversion. However, at intermediate 
rates the fluctuations are large, exceeding 50% of the mean number of bosons. We see from the 
graph that again the agreement between the exact and semiclassical curves is excellent and we 
expect it to be better for larger values of S. Moreover, one can show that the whole distribution 
function of the boson occupation number, which interpolates between the exponential distribution 



for fas t quenches and the Gumbe 



TWA (Altland et al 



2009; 



Kiegell . 



dist ribution at slow quenches, is accurately reproduced by 



2009) 



V. DERIVATION OF THE RESULTS WITHIN THE PATH INTEGRAL FORMALISM. 



A convenient framework connecting quantum and classical dynamics is given by the Feynman's 
path integral representation of the time evolution. Th en classical traj ectories appear in the form 



19941 ). As it is well known from 



of the saddle point of the action in the path integral (IShankar 
the standard quantum mechanics the classical limit is naturally recovered in the Heisenberg rep- 
resentation, where operators corresponding to various observables depend on time, while the wave 



funct ion (or the density matrix) only describes the initial state of the system (ILandau and Lifshita . 



1981 



). In this section we will show how combining path integrals and the Heisenberg represen- 



tation reproduces all the results of the previous sections. The formalism here will be somewhat 



resembling the one used in the derivat ion of the Keldysh diagrammatic technique (jKamenev 



2005 



Kamenev and Levchenko 



2002, 



20091 ) . However, there are important differences. Namely the 



Keldysh technique is usually derived in the Schrodinger representation and the elementary ob- 
jects there are the single-particle Green's functions. Also the majority of applications of the 
Keld ysh technique deal with s t eady state situations, where the initial conditions are unimpor- 



tant (Kamenev and Levchenko 



2009). While here we will work in the Heisenberg representation 



directly with the phase space variables and the initial conditions are of primary importance to 
us. The Heisenberg representation has another advantage over the Schrodinger representation that 
one does not have to calculate the density matrix, which is a tremendously complicated object in 
many-particle interacting systems containing exponentially large amount of information. Within 
the Heisenberg picture we aim to calculate expectation values of particular observables. 

Partially the results of Sec. IHII were obtained in literature earlier analyzing the von Neumann's 
equation for the density matrix in the Wigner form. We will briefly review this approach in the next 
section. The path integral formalism, in our opinion, gives a more natural and complete way to 
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obtain these results. Besides these earlier approaches failed to reproduce the structure of quantum 
corrections and non-equal time correlation functions. Another important goal of this section is to 
show that the Wigner function, Weyl ordering, Bopp operators and other related concepts naturally 
emerge in the path integral representation of time evolution without the need to make any extra 
assumptions or justify usage of these concepts a-posteriori. To keep the discussion relatively short 
we will outline only the key steps in the derivations skipping unessential intermediate details. 

A. Coordinate-momentum representation. 

1. Truncated Wigner Approximation. 



Let us start from a simple situation, where we are interested in the time evolution of the expec- 
tation value of some single-particle observable Q(p,x,t). According to general rules of quantum 
mechanics the latter is given by 



{(l(x,p,t)) = Tr \pT T e i/h tiK(T)dT^~^^ t ^ e -i/hftH(T)dT 



(161) 



whe re p is the (tim e independent) density matrix, T T stands for the Keldysh time order- 



ing (jKamenev 



2002), where later moments in time always appear closer to the observable SI. 



Next we split the evolution operator into the product: 

N 



r T e«Wr^ JJ e 



i/hH(n)AT 



(162) 



i=l 



where At = t/N and Tj = iAr. Finally we insert the resolution of identity \pi){pi\Xi){xi\ between 
all the terms in the product. At the end we will take the limit N —> oo. We use the standard 
convention for the overlaps between the coordinate and momentum basis states: 



(xi\pi) = exp[ipiXi/H], (pi\xi) 



1 



2irh 



exp[-ipiXi/h}. 



In this way we map the quantum evolution to the path integral over propagating classical phase 
space variables: x\ = x{t\) and pi = p{jj). Because we have two exponents in the path integral 
we have to introduce two sets of phase space variables: (x/i, Pfi) and (xi,i, Phi ), where indices 



f and b imply "forward" and " 



Kamenev and Levchenko 



2009; 



aackward". I t is int uitively clear (see also Refs. (jKamenev . 



Polkovnikov 



2002; 



2003al )) that the symmetric combination of the phase 



space variables: Xi = [xfi + xu[/2 and Pi = [pfi + pbi}/2 correspond to the classical coordinate and 
momentum, while their differences = xn — x^i and rji = pfi — pu describe quantum fluctuations. 
Indeed in the classical limit there is a unique trajectory corresponding to the fixed initial conditions, 
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while quantum mechanically one has to sum over all possible trajectories. Our goal is thus to 
develop a systematic expansion of (Cl(x,p,t)) in £j and rji, while treating classical variables Xi and 
Pi exactly. After inserting the complete basis and changing the variables to Xi,£i,pi,rji Eq. (|16ip 
assumes the form: 



{n(x,p,t)) 



dx dp Wo(po,xo) 



Dx{t)Dp{t)D^{t)Dt]{t) 



exp< 



dr 



., ,dp(r) . .dx(r) „, / . . r?(r) . . f(r) 



dr 

Uw \p(r) 



Or 



2 ' v ' 2 ' 

> n w (x(t),p(t),t) 



(163) 



where 1-Lw(p(t),x(t),t) is the Weyl symbol of the Hamiltonian 'H(x,p,t). The derivation of thi s 



expression is straightforward and very similar to the one described in Ref. (jPolkovnikov 



2003a|) 



in the coherent state basis. The expression in the brackets is nothing but the difference between 
forward ( Sf) and backw ard (Sb) classical actions appearing in the path integral of the evolution 



operator (IShankarl . 



1994 ): 



S 



f,b 



dr Xf tb (r 



dPf,b( T ) 



+ K(Xf,b(l~),Pf,b(T),T). 



(164) 



There are two subtle points in this derivation. First, one needs to carefully follow the bound- 
ary terms at r = and r = t. This can be straightforwardly done by working with the dis- 
crete time steps in the f unctional integral and taking the continuum limit at the very end (see 



Ref. ( Polkovnikov 



2003al ) for the analogous discussion in the coherent state basis). One can check 



that integrating over £(0) and r/(0) yields the Wigner transform of the density matrix (given by 
Eq. (0)), which is interpreted as the distribution function of the classical initial conditions xq and 
Pq. Integrating over £(t) and rj(t) results in the emergence of the Weyl symbol of the operator 
Q(x,p,t): ^lw( x (t),p(t),t), which is given by Eq. (JH). Remarkably both the Wigner function and 
the Weyl symbol appear automatically in the path integral. The second subtle point refers to 
the Weyl ordering of coordinates and momenta in the Hamiltonian in the exponent of Eq. (1163j) . 
For the Hamiltonians of the type of (|4ip which split into the sum of kinetic and potential en- 
ergies, this issue is not important because the Weyl ordering does not produce any additional 
terms. However, in a more general situation, where the Hamiltonian can contain cross-products 
like x 2 p 2 the issue of the correct ordering becomes relevant. Naively in this case one should use 
the "normal" ordering where x operators appear on the left of p operators because in this case 
{xi\exp[i'H(x,p,Ti)A.T]\pi) exp[i7i n (xi,pi + i, Tj)Ar], where n emphasizes that quantum-classical 
correspondence is taken in the Hamiltonian written in the normal form. However, this turns out 
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to be not accurate. The fact that the time indices in x and p in the Hamiltonian are not the same 
(or more precisely that the momentum p is evaluated at a slightly later time than the coordinate 
x) introduces the correction to the classical Hamiltonian % n modifying it to T-Cw We will discuss 
this issue in detail in the Appendix iBl 

Expression (I163P has in fact very transparent interpretation. The first integral is taken over 
all possible initial values of the classical coordinate and momentum xq and po weighted with 
the Wigner function. The latter thus plays the role of the probability distribution of the initial 
conditions. Next in Eq. (|163p there is a functional integral over classical and quantum fields 
with classical fields satisfying the initial conditions x(t = 0) = xq and p(t = 0) = po- And finally 
Q,w(x(t),p(t), t) is the Weyl symbol of the operator Cl(x,p, t), which expectation value is computed. 
As written, Eq. (|163p is formally exact. However, to make use of this representation one needs to 
find an efficient way of computing the functional integral over all trajectories. In this work we are 
focusing in expanding this integral in quantum fluctuations near the classical limit. In principle 
one can consider other expansions, e.g. in the interaction strength. The action in Eq. (|163p 
can also have other non-classical saddle points, which potentially describe such phenomena as 



tunneling (jKamenev and Levchenkd . 



20091 ). The latter can not be understood through the analytic 
expansion of the dynamics in powers of h. We already argued that 7/(t) and £(t) are responsible 
for the quantum fluctuations, so it is useful to expand the integrand in Eq. (I163p in powers of r/ 



and £. In the leading order we have 

r/(r 



H 



P{r) + 



U 



P{T) — ,X[T) — ,T 



^rr^ — + £( T ) 



(165) 



dp(r) dx(r) 
The derivatives of the Hamilt onian with respect to canon ical coordinate and momentum are just 



the classical Poisson brackets (Landau and Lifshitz 



1982 1: 



dnw(p,x,T) 

dx(r) 
dU w (p,x,T) 

dp(r) 



{U w (p,x,t),p} , 
-{H w {p,x,t),x} . 



(166) 



If we stop at this order of the expansion in r/(r) and £(t) then the functional integral in Eq. (|163p 
is trivial. Indeed, integrating out r/(r) and £(t) gives the 5- function constraint on p{r) and x(t): 



^ = {x,H w } , ^ = {p,7-Lw} 



(167) 



Thus taking the functional integral over the classical fields x(r) and p(j) becomes equivalent to 
simply solving the Hamiltonian equations of motion. Only the integral over the initial conditions 
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remains. So in the leading order in quantum fluctuations we recover the semiclassical (truncated 
Wigner) approximation given by Eq. (|44p . Straightforward generalization of this derivation to 
many-particle system immediately gives the same result with the only difference that now integra- 
tion is taken over coordinates and momenta of all particles in the system. 



2. Nonequal time correlation functions. 

Above we described how one can determine the expectation value of a single-time observable 
Q(x,p, t). We saw that one needs to use the Weyl symbol for the observable, which corresponds 
to first symmetrically ordering operators x and p and then substituting quantum operators by the 
corresponding classical variables. Now we will be interested in a more general case of finding the 
expectation value of the product of two (or more) observables evaluated at different times: 

{Cl 1 (x,p,t 1 )n 2 (x,p,t 2 )). (168) 

Let us first assume that t\ < t 2 . Later we will generalize our results to multiple time correlation 
functions and other types of ordering of the operat ors. In the c oherent state picture the non -equal 



2009 



I.Plimak et al. 



2001a! ) using 



time correlation functions were analyzed in Refs. (IBerg et al. 
related but somewhat different approach. 

Before proceeding with a formal derivation let us discuss where the potential complications come 
from. Classically the expectation value (I168p can be understood as average of the result of a joint 
measurement of first the quantity £li and then Q 2 . At the classical level (the ideal) measurement 
itself does not have any effect on the consequent dynamics thus one can write 

£ll{x,p,t\)VL 2 {x,p,t 2 ) = J J dxodpoW c i(x ,po)£li(x(ti),p(ti),ti)tt 2 (x(t 2 ),p(t 2 ),t 2 ) 

= J j duidu 2 uiu 2 P{ui,ti)P(u 2 ,t 2 \uji,ti), (169) 

where we use overline for the classical averaging to distinguish it from the quantum expectation 
value, 

P(ui,h) = J J <feo<^Wc{(xo,J%)ni(z(ti),p(ti),ti)%i-ni(a;(*i),p(ti),ti)) (170) 

is the probability that at the time t\ the measurement of J7i will give the value u\. Similarly 
P{uj 2 ,t 2 \uJi,ti) is the conditional probability that the measurement of £l 2 at the moment t 2 will 
give the value uj 2 given the measurement of f^i at t\ gave the value u\. Equation (|169p can be 
interpreted that the result of the measurement of £l\ simply restricts the possible initial conditions 
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xq and po to those which result in ili(x(ti),p(ti),ti) = uj\ but does not affect the consequent 
dynamics. In quantum mechanics the situation is more complicated. Indeed the very fact of 
measurement of ill has an effect on the consequent evolution and thus on the measurement of 
Mathematically it should be manifested in some form of back action at moment t\, which affects 
the consequent time evolution. The other important issue, which differentiates between quantum 
and classical multi-time measurements is the correct ordering of operators in the corresponding 
correlation functions. This discussion goes be yond the scope of this paper; it is available in li t eratu re 



for various specific situations (see e.g. Refs. 



Clerk et al 



2008; 



Nazarov and Kindermann 



20031 )1. 



We will just discuss how one can compute the correlation functions of the type fjl68f> within the 
phase space representation without connecting them to actual measurement procedures. 

According to Eq. @ the Weyl symbol corresponding to the product of two operators is formally 
defined as: 



(ni(ti)n 2 (t 2 )) 



w 



Air 



\ 



x(t x )+Z/2 



xn 2 (x,p,t 2 )e- l/h ^ H{r)d7 



-v/ 2 ) exp 



ni(x,p,t 1 )T T e i/h # H( - T) * r 



--[p{t l )x(t l ) + 



T]x(ti 



(171) 



It is easy to check that if Q 2 = 1 this expression indeed becomes equivalent to Eq. The 
equation (jlTlf) can be further analyzed similarly to Eq. (1161 j) . The result of this analysis is that 
at the moment t\ the phase space variables x(t) and p(t) become discontinuous acquiring jumps 
5x\ and 5p\ with the probability distribution, which depends on details of the operator Oi. If the 
latter is written in the ordered form where all x operators appear on the left of p operators the 
probability distribution of these jumps assumes the following form: 



Wn^xt^p!) 



d^drj 



(172) 



where we put the subindex Qi to highlight that these jumps depend on the form of the operator 
di. If Hi is the identity operator then the integral becomes a product of two 5- functions enforcing 
5x\ and Spi to be equal to zero. Then there is no jump in either x(t\) or p{t\) as it should be. If 
the operator il is written as a product of powers of x and p: il nm = x n p m , then the probability of 
the corresponding jump can be found using the generating function: 

d^drj 



R(a, f3, Sx, 5p) 



2?r 



-2ia(x+£/2)+2il3(p-ri/2) e 2i5p6x/h e i£6p/h+iri8x/h 



(173) 



To shorten notations we dropped explicit time indices for phase space variables x and p. Then 



Wn,m(Sx,8p) 



2 rH 



-d^dfR(a,P,Sx,Sp) 



(174) 
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The integral in Eq. (|173p can be easily evaluated and we find 

R{a, 0, Sx, Sp) = e -^x+2^ P +2i5x6p/h s ^ p _ a fy S ( Sx _ p h y ( 175 ) 

There are two particularly important cases n = l,m = corresponding to the coordinate operator 
(Oi o = x) and n = 0, m = 1 corresponding to the momentum operator (Oq i = p)- In the former 
case we find 

W 1)0 (6x, Sp) = e 2i5xS v' n (x - y 5(5p)5(5x). (176) 

Using standard rules of manipulating with (^-functions and noting that the derivative acting on the 
first exponent does not give any contribution we find 

W 1>0 (Sx,6p) = 6(5p)5(6x) (x + yj^) . (177) 

This expression is in turn equivalent to the representation (J46]): 

* (y ^ ( « + !a4r (178) 

From Eq. (|176p it is clear that the derivative here is understood as the response of the operator 
0,2 to the infinitesimal jump p{t\) — > p(t\) + 5p{t\). In Eq. (j46l) we used d/dp instead d/dSp to 
shorten notations. In the same way we find: 

W 0A (6x,6p) = 5(5p)S(5x) L-^JL} (179) 



leading to the substitution (pET|) : 



In the same spirit one can obtain the expression for the product xp: 

W ljl (Sx,Sp) = e 2 ^ ^ _ ( p _ s(Sp)S(Sx) (181) 

Now there is an additional subtlety that integrating 5-functions by parts we will get a nonzero 
contribution from the exponent coming from taking the double derivative with respect to 5x and 
5p so the result is (fT6|) : 

ih ih d ih d h 2 d 2 

xp^*p + ^- ^Wx + -J P W P + Tds^dSp- (182) 

In a similar fashion one can obtain representation of higher powers of x and p. However, there 
is a much simpler way to do this using just Eqs. ()178j) and (|180p . Note that in the derivation of the 
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representation of fii we never used any information about the second operator £l 2 . In particular, 
the same representation will be true if we have product of arbitrary number of operators ordered 
such that later time appears on the right e.g. ^1(41)02(^2)^3(^3) with t± < t 2 < £3 and so on. One 
simply needs to understand derivatives appearing in Eqs. (|178p . (I180p . and (11821) as acting on all 
operators appearing later in time (or on the right) of Then we can recover expression (|182l) 
from Eqs. (fTTH]) and (fT80|) noting that 

mm) = mm + ) - (*&) + § ^) (W + ) - y^y) ' (183) 

where ii+ stands for time approaching to t\ from above. Moving the derivative with respect to bp 
to the right and noting that we will pick up additional constant term ih/2 we recover Eq. (|182p . 
Similarly one can get the representation of other products of operators. E.g. p x = xp — ih. So the 
representation of px is given by Eq. (|182p with the negative sign in front of ih/2. The same result is 
also recovered by multiplying Eqs. (|180p and (|178|) in this order. We want to emphasize again that 
the same rule applies even in the case when Q, 2 is the identity operator. Then all derivative terms 
in Eqs. (JT75J), (|180p and (|182p drop and as a result we recover the Weyl symbol of the operator 
fil. This can be indeed checked by comparing the results of integration in Eq. Q and applying 
Eqs. (|178p and (|182p to specific operators. Thus for equal time correlation functions we just recover 
Bopp representation of the coordinate and the momentum which aut omatically reproduce Weyl 



symbols of arbitra ry operators as we discussed in Sec. HJ (see also Refs. (jBopp 



198 



Kubo . 



1961 



Hillery et all 



19641 )). It is notable that the Bopp representation automatically emerges from the 
path integral with the interpretation of the derivative terms as a response of future observables to 
the quantum jumps. 

In analyzing the correlation function (|168p we actually never used the fact that t% < t 2 : the time 
index in the path integral can increase or decrease since it is only the label for different phase space 
variables. So the correspondence (|46p and (|47p works if we are interested in correlation functions 
with the opposite ordering of times to (|168p . namely 

(Q 2 (x,p,t 2 )n 1 {x,p 1 t 1 )}. (184) 

It is convenient to assume again that t\ < t 2 but change the actual order of the operators so that 
the earlier one appears on the right. In this case using previous results one can reach the time 
t 2 and then simply propagate the equations of motion backwards in time until the moment t\. 
However, this representation is clearly not casual: quantum jumps at the later time t 2 affect the 
observable Q\ evaluated at an earlier time t\. It turns out that the casuality can be restored by 
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rewriting Eq. (|17ip in the following form 

d^drj I 



(n 2 (*2)fii(ti)) w 



4tt \ 



x(h) + £/2 



T JlKf t l HWtr^^ fe)e -*/»J$ H(r)d,T 



xOi(x,p,ti) 



p(ii) - r/ /2\ e - i P( tl ) :c ( tl )/ ri -^ p (* 1 )/( 2 ' l ) +i, ' a: (* 1 )/( 2S ). 



(185) 



This equation can be analyzed similarly to Eq. (117ip and the result is identical to (1172p with the 
only change that fei, 5p\ —> —Sxx, —8p\. In turn this leads to the change in sign in the derivatives 
in the representations (|178p and (|180j) . This change in sign is most easily interpreted as a result 
of simple calculus. Indeed the derivative with respect to 5x\ and 5p\ now acts on the operator 
f&2, which is on the left of Q.\. So Eqs. (j!78|) and ()180p have to be integrated by parts. We can 
combine the result for both types of orderings using the notations of left and right derivatives as 
in Eqs. (146p and (j47j) . The choice of the left or right derivative is arbitrary. However, it is unique 
if we want to preserve the casuality of our description. If earlier time appears on the left like in 
Eq. (|168p then we should use right derivatives in the representation of Vt\ and conversely in the 
ordering like in (1184p we have to use left derivatives. 

These results can be further generalized to multi-time correlation functions. At each moment of 
time we use one of the representations (|46p or (|47p depending on the convenience. As we saw in the 
case of the two-time correlation functions the choice of the representation is not unique. However, 
it becomes unique if we insist on the casuality. In general there are two special types of correlation 
functions (i) those which do not require any jumps so that one can use a simple quantum classical 
correspondence x — > x, p — > p and (ii) correlations which have casual representation. The first type 
of correlations can be directly calculated in TWA. The second type of correlations is very important, 
because only such correlations should appear in physical observables if we insist on casuality of the 
phase space representation of the quantum evolution, i.e. on that this representation is physical. 

The first type of correlations correspond to special, time s ymmetric, orde r ing of coordinate 



and momentum operators, which was first intro duced in R efs. 



somewhat different name and later used in Ref. (jBerg et al 



I.Plimak et al. 



2001a 



bj) under a 



20091 ). Let us first point that for two 



operators the time-symmetric ordering simply coincides with the symmetric ordering. Indeed if we 
are interested in the correlation functions of the type 



(x(ti)x(t 2 ) +x(t 2 )x(h)}, 



(186) 
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where we again assume t\ <t 2 , then according to Eq. ([36)) we should use the following substitution 
(x(h)x{t 2 ) + x(t 2 )x(h)) -> (x(h) + ^tt|-t1 x(t 2 ) + x(t 2 ) (x(h) ' h ° ^ 



a;(*i) + 



x(* 2 ) + 



2 d P {h) j 
ih d 



2 dp{ti 



x{t 2 ) = 2x(t 1 )x(t 2 ). 



(187) 



2 Spfa); 2 3p(^ 

Similarly one can check that the symmetric combination of x and p operators does not acquire any 
jumps: 



x(h)p(t 2 ) +p(t 2 )x(t 1 ) -> 2x(t 1 )p(t 1 



(188) 



This result can be straightforwardly generalized to larger number of operators noting that 

x(ti)A + A T x(ti) x(ti)(A + A T ), (189) 

where A is an arbitrary combination of x(tj) and p(tj), A T is the transposed operator where the 
same combination of operators x{tj) and p(tj) appears in the reverse order. E.g. if A = x{t 2 )p{t^) 
then A T = p{t%)x(t 2 ). Equation (|189p is still true if instead of x{t\) we will use p(t\). To avoid 
jumps appearing in the representation of A, the latter should have similar structure i.e. A = 
x{t 2 )B + B T x(t 2 ) or A = p(t 2 )B + B T p( t 2 ). The operator B s hould again have this form and 



so on. Following Refs. (IBerg et al 



2009; 



I.Plimak et al 



2001al ) we can define time-symmetric 



ordering symbol according to the following recursion relation: 



T[t(t 1 )t(t 2 ),..4(tn)} = \ U/,iT[i(/,) 



,£(tn)]+T[i(t2),...,i(tnMh) 



(190) 



2009J) time- 



where ^ stands for either x or p. We would like to emphasize that in Ref. (IBerg et al. 
symmetric ordering was defined with the additional requirement that t\ < t% < . . . , t n . As we will 
show next this enforces the casuality of the time evolution. Note that even though the operators x(t) 
and p{t) map to numbers x(t) and p(t) in the time-symmetrically ordered operators for arbitrary 
order of times ti, ... ,t n , one can not still reorder them i.e. the average of x(ti)x(t 2 )x(t3) is not 
equivalent to the average of x(t 2 )x(ti)x(t3). The reason is that the full representation of quantum 
evolution (beyond TWA) contains quantum jumps, which are sensitive to the ordering of x(tj). 
Only within the semiclassical (TWA) approximation are such permutations allowed because for 
each initial condition the classical trajectory is unique and one can arbitrarily reorder classical 
numbers. 

Now suppose that we are interested in a non-symmetrized correlation function of the type 



(191) 
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We want to understand for which orderings of time one can represent this correlator through 
quantum jumps in a casual way, i.e. in a way that jumps in x or p at moments t a only affect 
observables appearing at later moments in time. We saw above that if there are only two times 
involved then such representation is always possible. However, if we have a multi-time correlation 
function this is no longer true. Indeed let us consider a particular situation of a three time 
correlation function, e.g. 

(x(ti)x(t 2 )x(t 3 )) (192) 

Then for t\ < t 2 < t% there is an obvious casual representation where we use right derivatives in 
representation of both x{t\) and x(t 2 ): 

We note that in the leading order in h the jumps in the momentum at t± and t 2 are additive. In 
the next order in h 2 there is a double jump appearing in the form of the second order response of 
x(t 3 ) to p(ti) — > p{t\) + 5p(t\) and p(t 2 ) — > p(t 2 ) + Sp(t 2 ). Similarly for the opposite ordering of 
times t 3 < t 2 < t% one should use left derivatives resulting in 

where we reordered times in the casual way and removed arrows over the derivatives. 

There are two additional orderings for which there is a casual representation where one has to 
use combination of left and right derivatives. Thus if t\ < £3 < £2 we should use 

x{t x ) X {t 2 )x{t,) (x( h ) + j (x(t 3 ) - -— j x(t 2 ), (195) 

and finally for t 3 < t\ < t 2 the causal representation gives 

mnt 2 )x(t 3 ) (x( t3) - (x( tl) + f ^) x( t2) , (196) 

Note that the two remaining orderings where t 2 is the earliest time do not have any casual 
representation. Indeed whether we use the right derivative for x(t\) or left derivative for x(t 3 ) we 
will necessarily have to calculate the response of x(t 2 ) to the jump occurring at a later time. 

It is now straightforward to generalize this discussion to arbitrary multi-time correlation func- 
tions of type (I19ip . Casual representation exists only if the earlier time is never sandwiched 
between two later times. It is clear that there are 2 n orderings of times t%,t 2 , . . .t n satisfying this 
condition of casuality (compared to n\ of total number of orderings). Indeed to avoid sandwiching 
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the earliest time should appear either first or last in the sequence so there are only two distinct 
possibilities: either t\ or t n is the earliest time. Suppose that it is t\. Then there are only two 
possibilities for either t<i or t n being the next earliest time (if t n is the earliest time then either t\ 
or t n -\ should be the next earliest time) and so on. Interestingly only the orderings of correlation 



functions allowed by the casuality (so callec 



theory of measur ements 



Clerk et al. 



2008 



Schwinger-Keldysh orderings) ap pear in the quantum 



and (23) in Ref. (jNazarov and Kindermann , 



2003) (see e.g. Eqs. (22) 



Nazarov and KindermannL 

2OO3I )). This suggests that there arc potent hilly dee]) 



connections between casuality of the representation of the multi-time correlation function in the 
phase space and the possibility of measuring them. I.e. it is likely that only casual correlation 
functions are physical. 



3. Beyond truncated Wigner approximation. 

So far treating the exact expression (|163p for the evolution we considered only the leading linear 
terms in quantum fields £(r) and t](t) (see Eq. (|165p ). These terms resulted in the semiclassical 
(truncated Wigner) approximation. To get further quantum corrections one needs to include 
higher powers in the Taylor expansion of the Hamiltonian in quantum fields. For the quadratic 
(harmonic) Hamiltonians the higher order terms are absent and thus Eq. f|165j) and TWA are exact. 
For conventional non-relativistic Hamiltonians (141 p the kinetic energy is quadratic in p. Therefore 
higher order powers of quantum fields come only from the expansion in £. If one is interested in 
e.g. quantum corrections to the time evolution of a relativistic system, where the kinetic energy 



is of the form y m? + p 1 , then it is necessary to consider higher order terms in the field 77 in a 
similar fashion. Expanding Eq. (|165|) to the third order in £ and substituting this expansion into 
Eq. (fT63|) we find 

dx dpoWo(x ,po) J J Dx{t)D^{t) 

For the Hamiltonian (I4ip the momentum p(r) still satisfies the classical Hamiltonian equation of 
motion p(r) = mx(r) (because there are no cubic terms in rj) and thus the functional integral in 
Eq. (|197p is taken only over coordinates. This integral can not be evaluated exactly. However, 
one can find it perturbatively by expanding the exponent containing the third power of £ into the 
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exp (< fireugm + i /Vw^ + . . . {198) 



Taylor series: 



\24h Jo s v y 5x(r) 3 J 2AhJ s v y dx(^ 3 

Combining Eqs. ()197|) and (|198p and noting that 

1 d 3 1 5 

x 3 exp(ia:c) = -to-^t [5(a)] = — -rr8(a) 



2-7T i 3 3a 3 i 3 da 3 



we find 



/" /" / f l i h 3 d 3 V(x) d 3 \ 

(n((x),p,t)) « y J dx dp W (po,xo) [l - J o dr—-p Qx ^ 3 d6 ^ 3 j n w (x(t),p(t),t). 

(199) 

This expression is understood in the following sense. At some moment r the momentum p 
undergoes an infinitesimal quantum jump p{r) — > p(r) + 5p{r). Then the evolution proceeds 
according to the classical equations of motion. At the end one evaluates the nonlinear response 
of the observable to this quantum jump. Note that this jump indeed represents the quantum 
correction since it is accompanied by a factor of ti 2 . It is straightforward to generalize Eq. (|199p to 
include higher order terms in H. In particular, one can either have more quantum jumps occurring 
at different moments of time or one can have higher order quantum jumps corresponding to higher 
order terms in the Taylor expansion of Eq. (|165p . For example up to the fourth order in h one 
recovers Eq. ([52]) 

We can observe that formally Eq. (|199h can be rewritten through the quantum diffusion (see 
Eq. (|53p . which we repeat here for completeness): 



(Q(x,p,t)) « / / dxodp Wo(x ,po 

h 2 ^ d 3 v( 



n 



a- 



(x(t),p(t),i) + ^£^p) I d£F 3 (On w [x(t),p(t),t,8pi=tyZi\ \, (200) 



where Tj = it/N, At = t/N, Xi,pi = x(Ti),p(Ti). Formally we need to take the limit N — > oo or 
At — > 0. Equations (|199p and (I200p are equivalent if the function Fs(£), which can be interpreted 
as the probability distribution for the jumps, satisfies the following conditions: 

/oo /'OO rco rco 

F 3 ($di = o, / e^(CK = o, / e 2 ^ 3 (e)^ = o, / ^3(0^ = 1- ( 201 ) 
-OO J —CO J — CO J —CO 

The equivalence between the two representations of the correction (|199p and ()200p can be seen by 
expanding the observable £lw i n Eq. (|200p into Taylor series in 5pi\ 



9. 



w 



x(t),p(t),t,6pi = iyAn « n w [x{t),p(t),t\ H — (VAr 

J ddpi 



| l d 2 n w [x(t),p(t),t,5pj] 22/^-5 | l d 3 n w [x(t),p(t),t,5 Pi } ^^ | / Ar 4/3) / 202 n 
2 d5pf ^ 6 85p 3 
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Here all the derivatives are evaluated at dpi = 0. Now let us observe that the first three terms in 
this equation vanish upon integrating over £ weighted by F^{^) because the moments of £ up to the 
second are equal to zero (see Eq. (1201 [) ) . The terms containing powers of At higher than one vanish 
in the limit At — > 0. So the only term contributing to Eq. f|200|) from the Taylor expansion (|202p 
is the one containing the third order derivative with respect to 5pi. Combining this observation 
with the last of Eqs. (|20ip we see that Eq. (I200p is indeed equivalent to Eq. (|199p . Clearly the 
requirements (|20ip do not define the distribution function uniquely. One possible choice 

would be a Gaussian multiplied by a cubic polynomial (see Eq. (|57p ). Another possible choice for 
satisfying Eqs. (I20ip is the combination of (^-functions: 

m) = m-2) -*(g + 2) -2^ + 1) + 2^ + 1) (2Q3) 

This representation of F^(^) explicitly reduces Eq. (|200p to the finite difference representation of 
Eq. ([T991) . 

Similarly one finds that for the higher order jumps the distribution function must satisfy 



/oo 
f"F 2n+1 (£K = 0, m<2n, 
-oo 



oo 
oo 

f n+1 F 2n+1 (Z)dti = 1. (204) 



and the expansion takes the form of Eq. (|52p where higher order corrections in h either come from 
higher order jumps or from having several low order jumps. 

Multidimensional many-particle generalization. The expression of quantum evolution through 
response to infinitesimal jumps can be straightforwardly generalized to the situation where we deal 
with many degrees of freedom. For example instead of Eq. (|199p we should use Eq. (|58p . which we 
repeat here for completeness: 



(Cl(x,p,t)) « J J dx dp W {po,x ) 



* 1 h 2 d 3 V(x) d 3 . , 

*T-5i75!I -5Z *:L aZ TZ^Z 73 +- ^w(x(t), p(t), t). (205) 



o 



3!2 2 i 2 dxadx^dx^ dp a (r)dpp(T)dpj(T 
Here indices a,/3,j can run over different spatial and particle indices. For example for N particles 
in three dimensions they take the values x±, y±, Z\,X2, 2/2, %2-, ■ ■ ■ xn, Vn-> z n- This expression can be 
again rewritten through stochastic quantum jumps yielding Eq. (|59l) : 

(Ci(x,p,t)) « J j dx dp Wo(x ,p ) 



3 m cr(a,p,i) 



n w (x(t),p(t),t), (206) 
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where cr(a, ft, 7) stands for all inequivalent combinations of indices a, /?, 7 (e.g. all combinations 
with a < (3 < 7) and i ? a>j g )7 (^) is a function which should satisfy the conditions that all its moments 
smaller than two vanish and 

[ [Hdt m Z a ^F Q ^(£) = l. (207) 
A possible choice for the function F a ^ is given by Eqs. (foT1) - (l63l) : 

F m {i) = -J= (I - U\ e"^ 2 J] (208) 

Wfl = ^ = 1} ^ e-^/a U ^ m) , (209) 



^(0=^6-^^ n ( 2i °) 

Here the different labels a, /?, and 7 stand for different indices. 

The equivalence between representations of quantum corrections (I59p and (|58jl (or Eqs. (j206H 
and ([205 p ) can be established by Taylor expanding Eq. (|206p in powers of 5p Q (Tj) like in the 
single-particle case. A simple observation shows that only third order derivatives survive both 
the integration of £ Q and the limit At — > 0. Using that the functions F a p~ satisfy the property 
(|207p we recover that indeed Eq. (|206p reduces to Eq. (|205p in the limit At — > 0. The sum over 
non-equivalent permutations of indices a, f3, 7 in Eq. (|206p ensures that we avoid double counting 
of identical terms. 



B. Coherent state representation 

In this section we will show details of derivation of the results of Sec. (jlll.Bjl . The derivation 
closely mimics the one shown above in the number-phase representation. However, there are some 
distinct features as well. 



1. Truncated Wigner approximation 



First let us consider a single time expectation value of some operator fl(ip, ffi, t). To shorten 
notations we will skip the single-particle indices in the bosonic fields ip and ijr and reinsert them 
only when needed. The derivation essenti ally repeats t hat of the previous section. Partially the 



results were published previously in Ref. 



Polkovnikovl . 



2003al |. In the Heisenberg representation 
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this expectation value reads 



t)) = Tr [pT r e^o ?i(T)dT A(Vi,V» t ,t)e- i ^«W^ 



(211) 



Because in the coherent state picture the Planck's constant plays the mere role of conversion 
between time and energy units we set h = 1 throughout this section to simplify notations. Next 
we split the integral over time into a product as in Eq. (|162|) and insert the resolution of identity 
\tp(Ti))(jp(Ti)\ between different terms in this product. Then we get (cf. Eq. (|163p ): 



#0#$Wb(V>b,^o) / Dr/>(T)Dfl>*(T)DT,(T)DT}*(T) 



oxp<j I (It 



*/ , ^W(t) ( s n(r) , n*(r) 



rj* (r) 



JV(^(i),^*(t),i), 



(212) 



2 ' r v ' 2 

where J J dipdip* ■ ■ ■ = J J d?R.ip d^ip/it . . . Here Wo^ipo^tpo) is the Wigner transform of the initial 
density matrix (see Eq. (I26p ): 



drj Q dr]o (ip 



^0 + y } <' 



IV>ol 2 -jlwl 2 e |(??oV , o-J?oi/'5)_ 



(213) 



The quantities fiw(^>(t), and Hwiipir), ip*(r), r) are the Weyl symbol of the operator f2 

and the Hamiltonian "H (see Eqs. (|24p and (|25p). The emergence of the Weyl ordering of the 
Hamiltonian in Eq. (|213p is somewhat subtle (see related discussion in the coordinate-momentum 
representation after Eq. fjl63f) ) . We will discuss this issue in Appendix [B] 4 . In the leading order 
in quantum fluctuations we expand the integrand in Eq. (|212p up to the linear terms in rj. Then 
the functional integral over rj(t) enforces the J-function Gross-Pitaevskii constraint on the classical 
field ip(t): 

dH w (m,V(t),t) 



id t ip 



{m,H W (ip(t),r(t),t)} ( 



dip*(t) 

and we recover the semiclassical (truncated Wigner) approximation ([7T 



(214) 



2. Non-equal time correlation functions. 

The discussion for finding non-equal time correlation functions in the coherent state represen- 
tation closely mimics that of Sec. IV.A.2I in the coordinate-momentum representation. So we will 



We note that in Ref. (|Polkovnikovll2003al ) the issue of the Weyl rather than normal ordering of the Hamiltonian was 
overlooked. As we discussed in Sec. IIII.Bl see Eq. (|69p . for the Hamiltonians with two body uniform interactions 
this difference only affects the overall phase. 



77 



only give a brief r eview of the deriv ation. We note that the results presented in this section were 



discussed in Ref. (jBerg et al. 



2009) using somewhat different derivation. 



As in Sec. IV.A.2I we start from considering a two-time correlation function: 



(ni(^ T ,ti)JM^ T ,t 2 )> 



(215) 



Assuming that t\ < t<i and repeating the discussion of Sec. IV.A.2I we find that at the moment 
t\ the classical fields tp(ti) and ip*(ti) undergo quantum jumps: ip(ti) — > tp(ti) + ^(H) ~^ 
i/}*(ti) + Sip* described by the probability distribution which can be deduced from the generating 
function 

e d e* d 



R(5ip, 6i/>*, e, e*) = vr exp[e^* + e*ip - \e\ 2 /2] exp 



8{m^)8{%8^). (216) 



In this way £l nm = (ip^) n ip m produces the following probability distribution for the jumps 



(217) 



e=0 



As before the two most important cases are Clio = V* corresponding to 



and Oq l = ip corresponding to 



. + 1 d 
+-■ 



2d5^ 



Wb,iW,5^*) = 7r(^-i- ° 



2d8ip 



(218) 



(219) 



For the number density fin = we find 

1 ( , 5 9 



2 35^* 



1 



a 2 



<5(9ftty)<J(3<ty>). (220) 



55^ / 4 dSipdSip* 

As in the coordinate momentum case one can note that it is sufficient to establish the correspon- 
dence between operators and classical fields based on expressions for Wi o an d Wo \ and using the 
property of the (^-function: S'(x) = —5(x)d/dx: 



0t -». -0* 



1 9 



? , 1 5 
^ — > W + -- 



(221) 



2 (9(5-0 ' r ' r ' 2<9<5V>*' 
These are nothing but the coherent state Bopp operators ([33]) and (fM|) where the derivatives are 
understood as a response to the infinitesimal quantum jumps. If we reinsert the time label t% in 
the fields ip and if)* we will recover Eqs. (|72|) and (|73p . As we discussed earlier for finding a single 
time expectation value Bopp representation of ip and can be used to reproduce the Weyl symbol 
of an arbitrary operator. 
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The analysis of the correlation function with opposite time ordering, where the earlier time 
appears on the right, goes along the same lines as in Sec. IV.A.2I We thus will not repeat the 
derivation and only quote the final result that in this case one can use right derivatives with 
opposite signs as it is indicated in Eqs. ([72]) and (fT3|h This correspondence gives the casual phase 
space representation of the expectation value of ^2(^2)^1(^1) for t\ < t<i- 

Generalization of two-time correlation functions to multi-time correlation functions also goes 
completely along the lines of the discussion given in the previous section for the coordinate- 
momentum representation. In particular, the correspondence is very simple, if) — > tp, 'ijr — > tp*, for 
time-symmetrically ordered correlation functions given by Eq. (I190p . where now £ stands for either 
ip or ip* . For example: 



^(tiMWfo), (222) 



where [A, B] + stands for the anti-commutator of the operators A and B. Formally this substitution 
works for any ordering of times ti,t2,ts. However, only when t\ < ti < £3 the casuality of the 
description is preserved and one does not have to propagate equations of motion backwards in 
time. Generally (like in the coordinate momentum case) the casuality is preserved for operators 
appearing according to the Schwinger-Keldysh ordering, where the operator corresponding to an 
earlier time can never appear between the two operators evaluated at a later time. The casuality 
then corresponds to a particular choice of left or right derivatives at different times. For example 
for t\ < t<i < £3 we have 

- - ~) (m) - m), (223) 

while 

V(tMt*W 3 ) (Ml) - ig^j) (m) + igj^j) mi (224) 

C. Beyond the truncated Wigner approximation. 

Calculating quantum corrections to Eq. (|7ip also follows along the lines of the previous section. 
We need to expand the Hamiltonian in Eq. (|212p to higher orders in the quantum field r/. Note that 
for the nointeracting systems TWA is exact. For the Hamiltonians with two-body interaction (|64[) 
only the third power of rj appears in the expansion. To simplify the discussion we will focus on 
the derivation of quantum corrections in the system described by the Hubbard model with the 
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Hamiltonian (1661): 



Hhub = -Jj2^th + + E V A^i + 7 E U i3 4^Mi- ( 22 5) 



Then Eq. (|212|) assumes the following form: 



(iJ) 



expj i J dr 



di/j o di/joW (i/j o ,tl> 



0; 



Dip(T)Dip*(T)Dri(T)Dri*(T) 



-« — ^ h x I +Vj( r > 1 * — " 1 - 



<9r dipj(r) 



dr dijjj (t) 



+ 7 E ^ foiW<00 + ^(^)^(r)] |^-(t) 



(hi) 



n w (ip(t),r(t),t). 



(226) 



Next we expand the exponent containing cubic terms in r\j and to the leading order we find 



1 + I // r E U v [vUt)Mt) + Vi(r)rpt(r)} \ Vj (r 
.dipj(r) dHw 



Dip(T)Dip*(T)Dri(T)Dri*(T) 



x exp <i dr 



(id) 



<9r (r) 



( .^(r) , 



+ 



<9r dij)j (r) 



(227) 



nwr(^(t),^*(t),t). 



The functional integral over rj and rf can now be taken formally by introducing source terms into 
the path integral ([227]) : 



5 = exp 



E E M( T «H'( T <*) ~ Sj(T a )T]*(T a )) 

j T a =0 



(228) 



where to avoid dealing with functional derivatives we discretized the time: r a = aAr, a 
0, 1, . . . t/ At. As usual in the end we will take the limit At —> 0. Then we have 



dSj(T aJ 



-S, rjjira) 



ds*Ar a )' 



(229) 
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where all the derivatives are evaluated at Sj(r) = 0. Then we can rewrite Eq. (|227p in the following 



way: 



(A(W,t)> = I f </r/>„<ft/>SWoM>o,</>S) / / D*(t)B^*(t)Dij(t)Dij* 



■0i(T Q ) 



5 



V>*(r 



9 



dst(r a ) 



x exp < i At 

I r Q =0 



. ipj(r a + At) - ipjija) - Sj{r a ) 

2 \ 



^*(t q + At) - ty(r a ) - s*(r a ) 



At 



+ 



At 




(230) 



The source terms Sj(r) should be set to zero in the end of the calculation. Now one can integrate 
over rj(r) and t]*(t). As a result we get that the classical fields tpj(r) satisfy the Gross-Pitaevskii 
equations of motion with additional source terms, which play the role of infinitesimal shift of ipj{r): 
V'j(' r ) ~~ i J j{ T ) + s j( T )- Changing notations Sj(r a ) = 5ipj(T a ), taking the continuum limit At — > 0, 
and using the property of <5-function d s 5(s) = —5(s)d s we find that Eq. (|23ip reduces to 



i 



dif)od*l>oW (il)o,iPo) 
d 



a^-(r) 



9 



c) 2 



d5rl)j{T)d5ii)*AT) 



(231) 

nHr(^(t),^*(t),t), 



where at the time r the classical fields ipiij) and ipj(r) undergo infinitesimal jumps 5ipi(r) and 
5ipj(r) and later at the time t the nonlinear response of the observable f2 is evaluated with respect 
to these jumps. Note that Eq. (|23ip coincides with Eq. (|75p for a particular choice of the interaction 
Uijki = Uij5kjSn. Similarly higher quantum corrections to TWA appear in the form of multiple 
quantum jumps. 

As in the coordinate momentum case one can represent nonlinear response through stochastic 
quantum jumps. Then one arrives to Eq. (1 76 p . It is straightforward to check that indeed Eq. (|76D 
is equivalent to Eq. (I75P by Taylor expanding Eq. (176p in the powers of £j and noting that the 
terms up to the second order in this expansion give zero contribution to the integral because of the 
requirement (|77p . The fourth and higher order derivative terms disappear in the limit At — > 0. 
So the only terms in the Taylor expansion contributing to Eq. (I76p are the ones containing third 
order derivatives. These terms precisely reproduce Eq. ([73]) because of the condition ([77]) . As it is 
straightforward to check the functions ([78]) - ([ST]) give a possible choice for the distribution of jumps 
satisfying the conditions ([77]) . 
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VI. DERIVATION OF THE TRUNCATED WIGNER APPROXIMATION FROM THE VON 
NEUMANN'S EQUATION. 

The results of the previous section can be p artially recovered di rectly analyzing the von Neu- 
mann equation for the density matrix (see Ref. (jHillery et all . 11984 ) for details): 



ihp = [H, p\. 



(232) 



Taking the Weyl transform of both sides of the equation and using Eq. (|12p for the Moyal bracket 
we find that in the coordinate-momentum representation this equation becomes: 



W = {n w ,w} 



MB 



■ — fiw sin 
n 



W. 



(233) 



If we expand the Moyal bracket in the powers of h and stop in the leading order then the Moyal 
bracket reduces to the Poisson bracket and the von Neumann's equation (|233p reduces to the 
classical Liouville equation for the Wigner function: 



w « {n w ,w} 



(234) 



where the Wigner function plays the role of the classical probability distribution. The semiclassi- 
cal (truncated Wigner) approximation is recovered by noting that the classical trajectories are the 
characteristics of Eq. (|234p along which the Wigner function is conserved. Noting that the expec- 
tation value of any operat or is equal to the av erage of the corresponding Weyl symbol weighted 



19841 ) we indeed recover that Eq. (|234p is equivalent to 



with the Wigner function (jHillery et all 
Eq. (|44p . The interpretation of each trajectory in TWA, unlike in the purely classical limit, is not 



very straightforward for the reason that the Wigner function is not positive definite and can not be 
interpreted as a probability of a particular realization of the initial conditions. So the individual 
trajectory does not immediately represent a possible outcome of a single experiment corresponding 
to a particular realization of the initial conditions. 

The quantum corrections to Eq. (|234p appear because of the higher order term s in expansio n 
of the Moyal bracket in powers of H. Thus to the order of h 2 we find (see also Ref. (Zurekj, 120031 )) 



W 



n 2 



d 3 U 



dW 3 



2 2 3! dxjdxkdxi dpjdpkdpi 



+ 



(235) 



This partial differential equation is much harder to solve than the classical Liouville equation since 
the method of characteristics does not apply. One can attempt to solve it in the case of say a 
single degree of freedom, but then there is no advantage in this method over the direct solution of 
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the Schrodinger equation. Another possibility i s to w rite the evolution operator for this equation 
in the path integral form (see Ref. (|Berg et all . 120091 )) and try to treat this term perturbatively. 
In this case one should recover the results obtained in this paper. There are no other methods of 
efficiently simulating Eq. (|235p known in the literature. 

The truncated Wigner approximation also was independently developed in quantum optics as 



an approximate soluti on of the quantum L i ouvil 



coherent state basis ([Gardiner and Zoller . 



2004 



e equation for the Wigner function written in the 



Walls anc 



later was extended to the systems of interacting bosons 



many applications to cold atom systems (see e.g. Ref. (IBlakie et al 



Milburn 



Steel et al 



19941). This approximation 



1998) and recently found 



2003)). TWA immediately 



follows from the application of the coherent state Moyal brackets f|32[) to the von Neumann's 
equation for the density matrix (|232j) : 



ihW = {H W , W} MBC = 2U W sinh 



2 



W. 



(236) 



Expanding this Moyal bracket in the powers of A c (which is equivalent in the expansion in 1/N) 
to the leading order we find 



ihW m {H W ,W} C 



(237) 



recovering the classical Liouville equation for the density matrix. This equation can be solved using 
the method of characteristics which exactly coincide with the classical Gross-Pitaevskii equations: 

&Hw 



(238) 



These characteristics again represent trajectories along which the Wigner function is approximately 
conserved. 

If we go beyond TWA then we need to expand the coherent state Moyal bracket up to the third 
order in A c . Then we recover a third order partial differential equation analogous to Eq. (|235j) : 

d 3 U w d 3 W d 3 U w d 3 W 



(239) 



For the Hamiltonians with two-body interactions (or generally with not necessarily number con- 
serving interactions involving no more than quartic terms in the fields i/> and $') there are no 
higher than third order terms in the expansion of Eq. (12361) in powers of A c and thus Eq. 



is exact (see also Refs. (jBlakie et al. 



2008; 



Steel et al 



1998)). As in the coordinate- momentum 



representation this equation is hard to solve because the method of characteristics does not apply. 
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We note that the (probably unfortunate) name truncated Wigner approximation appeared in the 
literature for Eq. (I237D (which is also equivalent to Eq. (I7ip ) because at this level of approximation 
one truncates the Liouville equation (I239P to the linear order in derivatives with respect to ift and 
ip*. Thus this terminology is akin terming the Poisson bracket as the truncated Moyal bracket. 
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Appendix A: Implementation of quantum corrections for the sine-Gordon and Bose Hubbard models. 

In this section we give some details of how one can implement quantum corrections in practice. 
First let us consider the sine-Gordon model (see Sec. IIV.CI) . which is described by the Hamiltonian 
(|125p . At the semiclassical level one has to solve the equations of motion (|128j) subject to random 
initial conditions distributed according to the Wigner function (|126p . In the chosen setup the 
Wigner function is gaussian both at zero and at finite temperatures and there is no problem in 
sampling over the initial conditions. 

In order to implement the quantum correction in this case we used quantum jumps given by 
the second term in the brackets of Eq. (|53p proportional to h 2 . The function Fs(£) is given by 
Eq. (|57p . To implement this jump in practice at each step of the Monte-Carlo sampling over the 
initial conditions we also select a random time r and a random position j of the jump. Then at this 
time we evaluate V3J = — V(t)/3 3 sin(/30j(r)), randomly choose £ from a Gaussian distribution with 
zero mean and unit variance and the rest of the function i ? 3(£), namely l/2(£ 3 /3 — £), we multiply 
by Vsj and h? /2. And finally we transform 5rij — > 5rij + £j y/ At. Here At is a sufficiently small 
interval which can be adjusted to improve the efficiency of the simulation: larger At corresponds 
to faster convergence but worse accuracy. We checked that the results for At = 0.1 and At = 0.2 
are practically indistinguishable. 

Then the quantum correction term to the expectation value of say £1 = cos(/3^o(t)) is evaluated 
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as the average over many runs of the following quantity: 

tMF(r)/3 3 



8 



■ sin 



;(r))(r/3-£)co S (/^ (i))), 



(Al) 



where r is a random time uniformly distributed in [0, t), j is the random position of the jump: 
j € [0, M — 1), M is the spatial size of the system (M = 20 for our particular example), £ is the 
randomly distributed Gaussian variable with zero mean and unit variance. At the moment r we 
shift the density nj(r) — > nj(r) + ^At. At each Monte Carlo run we also randomly choose initial 
conditions distributed according to the Wigner function. 

Next let us describe implementation of the quantum corrections in terms of the nonlinear 
response, which we used in the problem of turning on interactions in the system of bosons in 
optical lattice (see Sec. IIV.D.3I) . In principle, the two implementations of corrections in terms of 
quantum jumps and nonlinear response are equivalent. We arbitrarily decided to use jumps for the 
sine-Gordon model and nonlinear response for the Hubbard model to illustrate that both methods 
work. Which one is more numerically efficient in which situation still remains to be investigated. 

As in the previous case TWA is obtained by solving the classical discrete Gross-Pitaveskii 
equations with the initial conditions distributed to the Gaussian distribution ()148p . The Weyl 
symbol for the Hamiltonian can be obtained e.g. by using Eqs. (fT3l) and (iT2j) : 



(bM*)l 2 -i) 2 



(A2) 



where the first sum is taken over all nearest neighbor pairs. 

To find the first quantum correction to the TWA result we use Eq. (175D . which in the case of 
Hubbard model takes the form 



<«(*)>! 



/ 



16 



oe\ de 2 



2 



+ 



o 2 



7iw(t,e 1 ,e 2 ) 



(A3) 



where ei and e 2 represent the real and imaginary parts of an infinitesimal jump of ipi at the moment 
t': 



(A4) 



Hw(t, ei , e 2 ) depends on e\ and e 2 through the fields ipj{t). To get Eq. (|A3P from Eq. ([75]) we used 
relations (|74|) because in practice it is easier to manipulate with real rather than comlex numbers. 
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Evaluation of the correction according to Eq. (jA3J) requires computing third order derivatives. This 
can be done using finite differences. E.g. for an arbitrary function f2(ei,e 2 ) we have 

<9 3 Q(ei, e 2 ) fi(2e 1} 0) - fi(-2e l5 0) - 2fl(e 1 ,0) + 2Q(-e 1 ,0) 



d 3 ft(ei, e 2 ) n(ei, e 2 ) + fi(ei, -e 2 ) - 0(-e x , e 2 ) - n(-ei, -e 2 ) - 2fi(ei, 0) + 2fi(- ei , 0) 



(A5) 



aeiSei ~ 2e ie | ' (A6) 

where ei and e 2 are small. As in the previous example the actual values of t\ and e 2 should be 
determined by the compromise: larger values give faster convergence but smaller accuracy. Similar 
expressions can be used to evaluate the remaining two derivatives. In practice one can evaluate 
this differences by time-evolving the original fields ipj as well as twelve additional perturbed fields 
which obtained by shifting ipi(t) in different ways: 

ipi,i(t) = ipi(t) + ei, tpi,2(t) = ipi(t) - ei, 

ipi,z(t) = M*) + 2e i> M*) = Mt) - 2ei, 
ipi,5(t) = 7pi(t) + ie 2 , ipi,e(t) = ipi(t) - ie 2 , 
Aj{t) = A(t) + 2ie 2 , 7p ita (t) = ipi{t) ~ 2ie 2 , 
A,9( f ) = ^*(*) + e i + ie 2, ^i,io(*) = A(t) - ei + ie 2 , 
A,n(t) = ipi(t) + ei - ie 2 , ^^(i) = ^(t) - ei - ie 2 , 

Then one can obtain desired finite differences using appropriately displaced solutions of the Gross- 
Pitaevskii equations (according to Eqs. (|A3j) . (IA5j) . and (jA6[) ) . As in the sine-Gordon case we 
can evaluate time integral and the sum by randomly choosing time if E [0, t) and position at each 
Monte-Carlo run and then multiply the result by the product of t and the number of sites. 

Appendix B: Emergence of the Weyl ordering of the Hamiltonian in the path integral representation 
of the evolution operator. 

The emergence of the Weyl form of the Hamiltonian both in Eq. (]163p and (I212p is 
some what subtle because usu ally path integrals are based on the normal ordering of opera- 



tors (jNegele and Orlandl . 



19881 ). In the coherent state picture the normal ordering corresponds 
to creation operators appearing on the left of annihilation operators. In the coordinate momentum 
representation by the normal ordering we will understand the one corresponding to the coordinate 
operators appearing on the left of the momentum operators. Note that usually the Hamiltonian is 
the sum of kinetic and potential energies the former being a function of momenta only and the latter 
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being the function of the coordinates. Then the issue of ordering never arises because normal and 
Weyl orderings are identical. This issue will become important if the Hamiltonian has cross terms 
like x 2 p which will emerge e.g. in the systems with nontrivial vector potential. In the coherent 
state picture the choice of the correct ordering can be very important. Even for the Hamiltonians 
with two-body interactions like (I64p the Weyl compared to normal ordering introduces additional 
terms (see Eq. (I68p ). As we argued in Sec. IIII.BI for spatially uniform two-body density-density 
interactions the Weyl ordering introduces an additional term proportional to the total number of 
particles, which is conserved in time and only affects an overall phase in the density matrix (and 
can be removed by a global gauge transformation). Such phase is usually unimportant unless we 
are interested in various interference phenomena between different systems. However, if the inter- 
action between particles is spatially dependent due to some external potential, the additional term 
in the Hamiltonian can not be eliminated by a simple gauge transformation and have to be taken 
into account. At the end of this Appendix we will illustrate this point with an explicit example. 

Let us discuss the emergence of the Weyl ordering in the coherent state representation. The 
discussion in the coordinate-momentum representation would be completely analogous. In the path 
integral representation of the Hamiltonian we have to deal with the evolution operator written as 
an infinite product (see Sec. IV.B]) : 

J]exp[±i7i(^ t ,V,r a )Ar]. (Bl) 

a 

For simplicity we suppressed all spatial indices. The two signs correspond to the forward and 
backward evolution. It is convenient to proceed with the Hamiltonian written in the normal- 
ordered form. Then by introducing the resolution of the identity as |V'(t q ))(V'(t q ) exp[— |^a| 2 ]| at 
each time step we end up with the products of the type: 



n 



a 



exp 



± Yl ^ T «) (V( T «+l) - ^( T ») + in n {^(T a ),lP(T a+1 ),T a )ATj 



(B2) 



Here the subindex n emphasizes that the substitution ^ — > ip* and ip — > ip is taken in the normal 
ordered Hamiltonian. The key observation we have to make is that the field ip* appears in the 
earlier time than the field tp. As we argued earlier, when we discussed non-equal time correlation 
functions, the path integral automatically produces Bopp operators for the fields ip and "ijp in 
any observable. The same situation happens in the Hamiltonian, which can be formally seen by 
treating the term exp [i% n (^*(r a ), ^(r a+ i), r a )Ar] at each time step as an observable. Then we 
have to use the substitution in the Hamiltonian, ffi—tt/)*— |c^. We remind that the derivative 
term comes from integrating over the quantum field 77 (see Sec. IV.Bp . for the details. As we 
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discussed in that section the derivative terms acting on the observables appearing at later times 
give raise to various quantum jumps. However, there is also an additional contribution from the 
jumps in ip on the Hamiltonian itself coming from the fact that in the normal ordered form tp 
appears at an infinitesimally later time than ip* (see Eq. (|B2|) ). As we discussed in Sec. IIII.Bl these 
derivative terms produce the Weyl symbol of the Hamiltonian. In other words the Weyl symbol 
of the Hamiltonian is produced by the self-action of the quantum jumps on the Hamiltonian. 
The same argument shows that the Weyl ordering should be used in the coordinate-momentum 
representation. 




FIG. 16 Energy (per site) dependence on the product 5t for a periodic Bose-Hubbard chain consisting of 
eight sites with the modulated interaction potential (see the text for details). The black solid curve is the 
result of the exact calculation. Red circles correspond to TWA with the Weyl form of the Hamiltonian used 
in the time evolution and and blue squares show the result of TWA where the normal ordered Hamiltonian 
is used instead. 

To illustrate the importance of the Weyl ordering we will consider an explicit example as in 
Sec. IIV.D.3I but with slight modification that the atoms are only interacting in off cites: 

« = -J^(^. +1+ ^t +i ^ ) + M £ (B3 ) 

j j=even 

The classical normal-ordered counterpart of this Hamiltonian appears by simply substituting — > 




FIG. 17 Same as in Fig. [16] except for S = 1. The black line is the exact result. The red circles and 
blue squares represent TWA and TWA + first quantum correction evaluated using the Weyl form of the 
Hamiltonian. 

ip* and ip — >• ip in the Hamiltonian above. The Weyl symbol for this Hamiltonian contains an 
additional term: 

Wwr = _j^(^ +1 + ^ i ) + M £ |^| 4 -C/(t) ^ |^f. (B4) 

j j'=even j'=even 

The extra term now does not commute with the Hamiltonian and gives raise to the physical 
difference in the time evolution. 

As in the example considered in Sec. IIV.D.3l we will start in the noninteracting case and assume 
that U(t) = Uq tanh<5i. We will use the same parameters as in Fig. [TOl i.e. eight sites, on average 
one particle per site, J = 1, Uq = 1 and 5 = 2.5. In Fig. [16] we show the comparison of the exact 
calculation and TWA with the Hamiltonian in the Weyl and normal-ordered forms. Clearly the 
Weyl representation (red circles) gives much better accuracy. In addition to a larger mistake using 
normal-ordered Hamiltonian also yields an unphysical result where the energy is not conserved 
even when the interaction almost stops changing in time. In Fig. [T7] we show that as in Fig. [10] 
the first quantum correction significantly improves the accuracy of TWA. We use a slightly smaller 
value of 5 = 1 in Fig. [17] compared to Fig. [16] to introduce a larger mistake in TWA so that the 
effect of the correction is more visible. 
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